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}
(136) 
_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 





1
(5) 
2
(1) 
3
(1) 
4
(3) 
5

6

7
(2) 
8
(2) 
9

10

11

12

13

14
(15) 
15
(11) 
16
(9) 
17
(1) 
18

19

20

21
(1) 
22
(9) 
23
(5) 
24
(1) 
25
(3) 
26
(5) 
27
(5) 
28
(21) 
29
(9) 
30
(3) 

From: Roy Stogner <roystgnr@ic...>  20100428 23:13:13

On Wed, 28 Apr 2010, liang cheng wrote: > Stack frames: 10 > 0: print_trace(std::ostream&) > 1: QBase::n_points() const > 2: assemble_ldm(EquationSystems&, std::string const&) > > This is the content in the traceout_0_?????.txt, the ????? varies according to > different case. I am not good at debugging, the falling one seems the No.10, > which is the name of executable file '/dld3d_v5dbg() [0x441d69] '. Does that mean > the whole program collapsed because of something wrong in the code? No, the failure was the n_points() call in assemble_ldm(). The quadrature rule in that function isn't getting initialized properly. > /usr/include/c++/4.4/debug/vector:272:error: attempt to subscript container > with outofbounds index 0, but container only holds 0 elements. > > Objects involved in the operation: > sequence "this" @ 0x0x7fffd1904378 { > type = NSt7__debug6vectorIdSaIdEEE; > } > > No traceout_0_xxx files generated this time, I guess the error comes from the > vector definition, I will check the whole vectors later. This is a libstdc++ generated error; to get a stack trace there the best way is to run from within gdb and type "where" after it crashes. > The model is about the solid elasticity problem, which is a second order > time dependent PDEs, so probably the ex8 is a good example for me. Yes, I think so. > maybe I should get the steady case work before jumping to the transient problem. Probably a good idea. > > Yes. With one variable on a grid of a million trilinear hexahedra, > > for instance, the SparseMatrix will be a (roughly) 1000000x1000000 > > matrix (which fits into memory because we only store the 1000000x27 > > potentially nonzero entries) while the DenseMatrix will be a 8x8 > > matrix with all 64 entries typically nonzero. > > Sounds the Densematrix is quit efficient, Yes, for one element. A dense matrix for the whole problem would have 1000000x1000000 entries and would require around 8,000GB of memory. > I am thinking how difficult to develop another 2nd time rate solver > such as HHT. is it possible to get new time solver done with one > source file and header file? Yes, but difficult. I wouldn't attempt it yet. Using ex8 as a template for your own application isn't exactly easy, but it's a lot easier than using newmark_system.* as a template for your own integrator would be.  Roy 
From: liang cheng <goeasyon@gm...>  20100428 22:50:50

> You're using a quadrature rule somewhere, though: it's QBase::n_points() that's failing. In ex8 there are two calls to > n_points()  did you add a third, and in any case which of them is > failing? (look for the traceout files with stack traces, or run in a > debugger) For the failing one, how did you declare the rule? Is it > being initialized before use? Stack frames: 10 0: print_trace(std::ostream&) 1: QBase::n_points() const 2: assemble_ldm(EquationSystems&, std::string const&) 3: System::user_assembly() 4: System::assemble() 5: ImplicitSystem::assemble() 6: NewmarkSystem::assemble() 7: main 8: __libc_start_main 9: ./dld3d_v5dbg() [0x441d69] This is the content in the traceout_0_?????.txt, the ????? varies according to different case. I am not good at debugging, the falling one seems the No.10, which is the name of executable file '/dld3d_v5dbg() [0x441d69] '. Does that mean the whole program collapsed because of something wrong in the code? Yes, they do. These are the global sparse matrices and forcing vector > that define the ex8 algebraic problem on the entire domain. But they > are constructed from components like: > > > DenseSubMatrix<Number> >> Kuu(Ke), Kuv(Ke), Kuw(Ke), >> Kvu(Ke), Kvv(Ke), Kvw(Ke), >> Kwu(Ke), Kwv(Ke), Kww(Ke); >> > > which are the dense matrices (and submatrices, for each variable) on > each element. In ex8, stiffness.add_matrix(Ke, dof_indices) adds a > permutation (from local to global indices) of the Ke matrix to the > stiffness matrix. > > I didn't know this method, thanks for telling me this to convey the [Ke] to the [stiffness]. When I add these commands to the code, the previous error disappeared and another shows like /usr/include/c++/4.4/debug/vector:272:error: attempt to subscript container with outofbounds index 0, but container only holds 0 elements. Objects involved in the operation: sequence "this" @ 0x0x7fffd1904378 { type = NSt7__debug6vectorIdSaIdEEE; } No traceout_0_xxx files generated this time, I guess the error comes from the vector definition, I will check the whole vectors later. > > Not just required for the NewmarkSystem solver, but for the underlying > time integration algorithm... if you didn't want to use the Newmark > method then ex8 might not have been the best choice of starting point > for you. > > > The model is about the solid elasticity problem, which is a second order time dependent PDEs, so probably the ex8 is a good example for me. maybe I should get the steady case work before jumping to the transient problem. > > Yes. With one variable on a grid of a million trilinear hexahedra, > for instance, the SparseMatrix will be a (roughly) 1000000x1000000 > matrix (which fits into memory because we only store the 1000000x27 > potentially nonzero entries) while the DenseMatrix will be a 8x8 > matrix with all 64 entries typically nonzero. >  > Sounds the Densematrix is quit efficient, so why not the Newmark solver use the Densematrix also, Is it limited by the PETSc solver? I am thinking how difficult to develop another 2nd time rate solver such as HHT. is it possible to get new time solver done with one source file and header file? I greatly appreciate your time and clear explanation, it is very helpful. Thanks! Liang 
From: Roy Stogner <roystgnr@ic...>  20100428 21:48:08

On Wed, 28 Apr 2010, liang cheng wrote: > No quadrature rule was added in the fill_dirichlet_bc() because no such commands in the > ex8 neither, I don't dare to violate the example too much, so didn't plug the quadrature > rule in fill_dirichlet_bc(). You're using a quadrature rule somewhere, though: it's QBase::n_points() that's failing. In ex8 there are two calls to n_points()  did you add a third, and in any case which of them is failing? (look for the traceout files with stack traces, or run in a debugger) For the failing one, how did you declare the rule? Is it being initialized before use? > I guess this problem is related to the matrix definition, Almost certainly not, but let's try and clear up that too. ;) > maybe my guessing is incorrect. > The matrix system was defined as followed, > > SparseMatrix<Number>& stiffness = large_deform_system.get_matrix("stiffness"); > SparseMatrix<Number>& damping = large_deform_system.get_matrix("damping"); > SparseMatrix<Number>& mass = large_deform_system.get_matrix("mass"); > NumericVector<Number>& force = large_deform_system.get_vector("force"); > > However these matrix doesn't participate in computation. Yes, they do. These are the global sparse matrices and forcing vector that define the ex8 algebraic problem on the entire domain. But they are constructed from components like: > DenseSubMatrix<Number> > Kuu(Ke), Kuv(Ke), Kuw(Ke), > Kvu(Ke), Kvv(Ke), Kvw(Ke), > Kwu(Ke), Kwv(Ke), Kww(Ke); which are the dense matrices (and submatrices, for each variable) on each element. In ex8, stiffness.add_matrix(Ke, dof_indices) adds a permutation (from local to global indices) of the Ke matrix to the stiffness matrix. > These matrix seems required for the Newmark solver as I dig into the > newmark_system.C/h. Not just required for the NewmarkSystem solver, but for the underlying time integration algorithm... if you didn't want to use the Newmark method then ex8 might not have been the best choice of starting point for you. > Is this 'SparseMatrix' data structure quit different with the that > of 'DenseMatrix<Number>'? Yes. With one variable on a grid of a million trilinear hexahedra, for instance, the SparseMatrix will be a (roughly) 1000000x1000000 matrix (which fits into memory because we only store the 1000000x27 potentially nonzero entries) while the DenseMatrix will be a 8x8 matrix with all 64 entries typically nonzero.  Roy 
From: liang cheng <goeasyon@gm...>  20100428 20:09:14

> > > And what about the boundary fe object? Did you also reinit that and > with the correct side? > >  > John > Yes, I did the BC like the ex8, the Neumann BC is defined in the assemble() and the Dirichlet BC is defined in fill_dirichlet_bc(). Is it critical if no fe_face>attach_quadrature_rule() and reinit() in the fill_dirichlet_bc()? No quadrature rule was added in the fill_dirichlet_bc() because no such commands in the ex8 neither, I don't dare to violate the example too much, so didn't plug the quadrature rule in fill_dirichlet_bc(). I guess this problem is related to the matrix definition, maybe my guessing is incorrect. The matrix system was defined as followed, SparseMatrix<Number>& stiffness = large_deform_system.get_matrix("stiffness"); SparseMatrix<Number>& damping = large_deform_system.get_matrix("damping"); SparseMatrix<Number>& mass = large_deform_system.get_matrix("mass"); NumericVector<Number>& force = large_deform_system.get_vector("force"); However these matrix doesn't participate in computation. These matrix seems required for the Newmark solver as I dig into the newmark_system.C/h. Is this 'SparseMatrix' data structure quit different with the that of 'DenseMatrix<Number>'? I found the 'DenseSubMatrix' could be defined while the 'SparseSubMatrix' not. Since the governing equation here require the matrix structure like DenseSubMatrix<Number> Kuu(Ke), Kuv(Ke), Kuw(Ke), Kvu(Ke), Kvv(Ke), Kvw(Ke), Kwu(Ke), Kwv(Ke), Kww(Ke); So the question is building Sub Matrix for the SparseMatrix like 'stiffness '. Maybe I mix the two separated problems. I have double checked the codes with the ex8, the only difference up to now is the submatrix system. For the variable of wave equation in ex8 is only one (pressure), while here in my case is three variables, which leads to a 3x3 stiffness matrix. I thought the solid element is easy, but it turns out almost kill me. I appreciate your time and answer. Thanks! Liang ******************************************************************** { // boundary condition below for (unsigned int s=0; s<elem>n_sides(); s++) if (elem>neighbor(s) == NULL) { * AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); QGauss qface (dim1, fe_type.default_quadrature_order()); fe_face> attach_quadrature_rule (&qface);* const std::vector<std::vector<Real> >& phi_face = fe_face>get_phi(); const std::vector<Real>& JxW_face = fe_face>get_JxW(); * fe_face>reinit(elem,s);* for (unsigned int qp=0; qp<qface.n_points(); qp++) { // Neumman boundary condition for (unsigned int i=0; i<phi_face.size(); i++){ Fu(i) += JxW_face[qp] * q_n * phi_face[i][qp]; Fv(i) += JxW_face[qp] * q_n * phi_face[i][qp]; Fw(i) += JxW_face[qp] * q_n * phi_face[i][qp]; } } // end face node loop } // end if (elem>neighbor(side) == NULL) } // end BC /******************************************************************* ******************************************************************** *** Apply the dirichlet boundary condition u, v , w *** ******************************************************************** *******************************************************************/ void fill_dirichlet_bc(EquationSystems& es, const std::string& system_name) { const Real y_coo = 6.; //BC coordinates libmesh_assert (system_name == "Disp_finite_deform_3d"); NewmarkSystem & system = es.get_system<NewmarkSystem> (system_name); SparseMatrix<Number>& matrix = *system.matrix; NumericVector<Number>& rhs = *system.rhs; const MeshBase& mesh = es.get_mesh(); const bool do_for_matrix = es.parameters.get<bool>("Newmark set BC for Matrix"); const Real current_time = es.parameters.get<Real>("time"); unsigned int n_nodes = mesh.n_nodes(); for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++) { const Node& curr_node = mesh.node(n_cnt); const Real penalty = 1.e10; // The penalty parameter. if (fabs(curr_node(2)  0.0) < TOLERANCE) { unsigned int dn = curr_node.dof_number(0,0,0); Real u_value = .0; Real v_value = .0; Real w_value = .0; rhs.add(dn, u_value*penalty); rhs.add(dn, v_value*penalty); rhs.add(dn, w_value*penalty); if (do_for_matrix) matrix.add(dn, dn, penalty); } // force or displacement if (fabs(curr_node(2)y_coo) < TOLERANCE) { unsigned int dn = curr_node.dof_number(0,0,0); Real w_value; if (current_time < .002 ) w_value = 0.1; // sin(2*pi*current_time/.002); else w_value = 0.2; rhs.add(dn, w_value*penalty); if (do_for_matrix) matrix.add(dn, dn, penalty); } } // loop n_cnt } 
From: John Peterson <peterson@cf...>  20100428 19:27:13

On Wed, Apr 28, 2010 at 1:50 PM, liang cheng <goeasyon@...> wrote: > On Wed, Apr 28, 2010 at 12:49 PM, Roy Stogner <roystgnr@...>wrote: > >> >> On Wed, 28 Apr 2010, liang cheng wrote: >> >> I made a 3D finite deformation elasticity input code and got the >>> executable >>> file, while the logical error displayed >>> in the screen when running it. >>> >>> Assertion `!_points.empty()' failed. >>> [0] /home/liang/Desktop/libmesh/include/quadrature/quadrature.h, line 116, >>> >> >> This usually means that the quadrature rule wasn't properly >> reinitialized on an element. Did you attach it to your FE object, and >> did you reinit() the latter? And what about the boundary fe object? Did you also reinit that and with the correct side?  John 
From: liang cheng <goeasyon@gm...>  20100428 18:50:37

On Wed, Apr 28, 2010 at 12:49 PM, Roy Stogner <roystgnr@...>wrote: > > On Wed, 28 Apr 2010, liang cheng wrote: > > I made a 3D finite deformation elasticity input code and got the >> executable >> file, while the logical error displayed >> in the screen when running it. >> >> Assertion `!_points.empty()' failed. >> [0] /home/liang/Desktop/libmesh/include/quadrature/quadrature.h, line 116, >> > > This usually means that the quadrature rule wasn't properly > reinitialized on an element. Did you attach it to your FE object, and > did you reinit() the latter? >  > Roy > Yes, I did it in the assemble subroutine as below. The system.init() is added in the man(). I am using both the [SparseMatrix<Number>& stiffness] and [DenseMatrix<Number> Ke] to define the stiffness matrix, actually only [DenseMatrix<Number> Ke] is used for computation, and the compiler warned me that the [SparseMatrix<Number>& stiffness] is not used, but still pass the compile. I am thinking if the Newmark solver requires the [SparseMatrix<Number>& stiffness] to define the [K] matrix instead of [DenseMatrix<Number> Ke], and the same thing for the [C],[M] matrix. Is this possible reason for the error? Thanks! Liang ************************************* void assemble_ldm (EquationSystems& es, const std::string& system_name) { libmesh_assert (system_name == "Disp_finite_deform_3d"); const MeshBase& mesh = es.get_mesh(); const unsigned int dim = mesh.mesh_dimension(); const Real G_neoHoo = es.parameters.get<Real>("shear modulus"); const Real K_neoHoo = es.parameters.get<Real>("bulk modulus"); const Real mass_dens = es.parameters.get<Real>("mass density"); NewmarkSystem & large_deform_system = es.get_system< NewmarkSystem> (system_name); const unsigned int u_var = large_deform_system.variable_number ("u"); const unsigned int v_var = large_deform_system.variable_number ("v"); const unsigned int w_var = large_deform_system.variable_number ("w"); FEType fe_type = large_deform_system.variable_type(u_var); // n1_var=0 SparseMatrix<Number>& stiffness = large_deform_system.get_matrix("stiffness"); SparseMatrix<Number>& damping = large_deform_system.get_matrix("damping"); SparseMatrix<Number>& mass = large_deform_system.get_matrix("mass"); NumericVector<Number>& force = large_deform_system.get_vector("force"); SparseMatrix<Number>& matrix = *large_deform_system.matrix; DenseMatrix<Number> zero_matrix; AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); QGauss qrule (dim, fe_type.default_quadrature_order()); QGauss qface (dim1, fe_type.default_quadrature_order()); * fe> attach_quadrature_rule (&qrule);* fe_face> attach_quadrature_rule (&qface); const std::vector<Real>& JxW = fe>get_JxW(); const std::vector<std::vector<Real> >& phi = fe>get_phi(); const std::vector<std::vector<RealGradient> >& dphi = fe>get_dphi(); const DofMap& dof_map = large_deform_system.get_dof_map(); DenseMatrix<Number> Ke, Ce, Me; DenseVector<Number> Fe; /********* 3 variables => 3x3 K matrix and 3 vectors *********/ DenseSubMatrix<Number> Kuu(Ke), Kuv(Ke), Kuw(Ke), Kvu(Ke), Kvv(Ke), Kvw(Ke), Kwu(Ke), Kwv(Ke), Kww(Ke); DenseSubMatrix<Number> Cuu(Ce), Cuv(Ce), Cuw(Ce), Cvu(Ce), Cvv(Ce), Cvw(Ce), Cwu(Ce), Cwv(Ce), Cww(Ce); DenseSubMatrix<Number> Muu(Me), Muv(Me), Muw(Me), Mvu(Me), Mvv(Me), Mvw(Me), Mwu(Me), Mwv(Me), Mww(Me); DenseSubVector<Number> Fu(Fe), Fv(Fe), Fw(Fe); . . ......define additional vectors segment . . . for ( ; el != end_el; ++el) { const Elem* elem = *el; dof_map.dof_indices (elem, dof_indices); dof_map.dof_indices (elem, dof_indices_u, v_var); dof_map.dof_indices (elem, dof_indices_v, u_var); dof_map.dof_indices (elem, dof_indices_w, w_var); const unsigned int n_dofs = dof_indices.size(); const unsigned int n_u_dofs = dof_indices_u.size(); const unsigned int n_v_dofs = dof_indices_v.size(); const unsigned int n_w_dofs = dof_indices_w.size(); * fe>reinit (elem);* Ke.resize (n_dofs, n_dofs); Ce.resize (n_dofs, n_dofs); Me.resize (n_dofs, n_dofs); zero_matrix.resize (n_dofs, n_dofs); Fe.resize (n_dofs); 
From: Liang <goeasyon@gm...>  20100428 18:21:58

> > Yes, if you use a separate ExplicitSystem. You don't want to add > variables to your main System unless you're planning to solve for > them. >  > Roy > Thanks Roy! Liang 
From: Boyce Griffith <griffith@ci...>  20100428 17:47:23

On 4/28/10 1:38 PM, Roy Stogner wrote: > > On Wed, 28 Apr 2010, Boyce Griffith wrote: > >> On 4/28/10 1:26 PM, Roy Stogner wrote: >>> >>> On Wed, 28 Apr 2010, Boyce Griffith wrote: >>> >>>> For AMR meshes, should these sorts of constraints be applied only for >>>> the nodes of level 0 elements? Or is this something that I don't have >>>> to worry about explicitly? >>> >>> You need to apply your constraints for all dofs that aren't already >>> constrained  everything except hanging nodes in adaptive h refinement >>> and their equivalent highp dofs in adaptive p refinement. In your >>> case it shouldn't matter whether or not you constrain hanging nodes as >>> well, but if you do then turn off forbid_constraint_overwrite when >>> adding them; by default we assume adding two different constraints to >>> the same node is a bug. >> >> Looking ahead to inhomogenous boundary conditions  > > A patch for that would be very appreciated, even if it only worked for > LAGRANGE variables at first, as long as the API was general enough to > be extensible later. This has *forever* been on our list of ideas > that are important/interesting enough to design but not quite > important/urgent enough to implement. It seems like it should suffice to allow a nonzero righthandside for the constraint equation. Am I missing anything? >> is there a mechanism for determining if a particular node is a hanging >> node, so that I can avoid reconstraining those DOFs? > > DofMap::is_constrained_dof() Each time that user_constrain() is called, have the constraints already been reset? I.e., in user_constrain(), is it safe to assume that the only constraints are hanging nodes (or their prefinement equivalent) and periodic boundaries?  Boyce 
From: Roy Stogner <roystgnr@ic...>  20100428 17:38:46

On Wed, 28 Apr 2010, Boyce Griffith wrote: > On 4/28/10 1:26 PM, Roy Stogner wrote: >> >> On Wed, 28 Apr 2010, Boyce Griffith wrote: >> >>> For AMR meshes, should these sorts of constraints be applied only for >>> the nodes of level 0 elements? Or is this something that I don't have >>> to worry about explicitly? >> >> You need to apply your constraints for all dofs that aren't already >> constrained  everything except hanging nodes in adaptive h refinement >> and their equivalent highp dofs in adaptive p refinement. In your >> case it shouldn't matter whether or not you constrain hanging nodes as >> well, but if you do then turn off forbid_constraint_overwrite when >> adding them; by default we assume adding two different constraints to >> the same node is a bug. > > Looking ahead to inhomogenous boundary conditions  A patch for that would be very appreciated, even if it only worked for LAGRANGE variables at first, as long as the API was general enough to be extensible later. This has *forever* been on our list of ideas that are important/interesting enough to design but not quite important/urgent enough to implement. > is there a mechanism > for determining if a particular node is a hanging node, so that I can avoid > reconstraining those DOFs? DofMap::is_constrained_dof()  Roy 
From: Boyce Griffith <griffith@ci...>  20100428 17:34:20

On 4/28/10 1:26 PM, Roy Stogner wrote: > > On Wed, 28 Apr 2010, Boyce Griffith wrote: > >> For AMR meshes, should these sorts of constraints be applied only for >> the nodes of level 0 elements? Or is this something that I don't have >> to worry about explicitly? > > You need to apply your constraints for all dofs that aren't already > constrained  everything except hanging nodes in adaptive h refinement > and their equivalent highp dofs in adaptive p refinement. In your > case it shouldn't matter whether or not you constrain hanging nodes as > well, but if you do then turn off forbid_constraint_overwrite when > adding them; by default we assume adding two different constraints to > the same node is a bug. Looking ahead to inhomogenous boundary conditions  is there a mechanism for determining if a particular node is a hanging node, so that I can avoid reconstraining those DOFs?  Boyce 
From: Roy Stogner <roystgnr@ic...>  20100428 17:26:59

On Wed, 28 Apr 2010, Boyce Griffith wrote: > For AMR meshes, should these sorts of constraints be applied only for the > nodes of level 0 elements? Or is this something that I don't have to worry > about explicitly? You need to apply your constraints for all dofs that aren't already constrained  everything except hanging nodes in adaptive h refinement and their equivalent highp dofs in adaptive p refinement. In your case it shouldn't matter whether or not you constrain hanging nodes as well, but if you do then turn off forbid_constraint_overwrite when adding them; by default we assume adding two different constraints to the same node is a bug.  Roy 
From: Boyce Griffith <griffith@ci...>  20100428 17:19:01

On 4/28/10 1:04 PM, Roy Stogner wrote: > > On Tue, 27 Apr 2010, Boyce Griffith wrote: > >> After adding DofConstraintRow's to the dof map, what generally needs >> to be (re)initialized? Is it a good idea to call >> EquationSystems::reinit()? > > The difficulty is in how you add the rows  it needs to be done in > user_constrain(), which gets called by EquationSystems::reinit(). If > you do it any earlier your rows are probably going to be deleted when > we construct the geometric/periodic constraints, and if you do it any > later then your rows might not be taken into account properly during > expansion of constraints. Aha  I missed that function. For AMR meshes, should these sorts of constraints be applied only for the nodes of level 0 elements? Or is this something that I don't have to worry about explicitly? Thanks,  Boyce 
From: Roy Stogner <roystgnr@ic...>  20100428 17:06:37

On Mon, 26 Apr 2010, Liang wrote: > Thanks! I thought the System::add_vector() stores the arbitrary vectors. > I appreciate your explanation. So the System::add_variable() will be a > good choice to define the flux component such as gradT_x, gradT_y, > If I just want to dump the these flux solution, right? Yes, if you use a separate ExplicitSystem. You don't want to add variables to your main System unless you're planning to solve for them.  Roy 
From: Roy Stogner <roystgnr@ic...>  20100428 17:05:04

On Tue, 27 Apr 2010, Boyce Griffith wrote: > On 4/27/10 10:20 PM, Roy Stogner wrote: >> >> On Tue, 27 Apr 2010, Boyce Griffith wrote: >> >>> I am trying to use the DofMap to impose Dirichlettype boundary >>> conditions within an explicit nonlinear mechanics solver which is >>> embedded in a larger fluidstructure interaction solver. I constrain >>> the velocity and force DOFs at Dirichlet boundaries to be zero by adding >>> a DofConstraintRow for each of the constrained DOFs. >> >> This should be a good method for homogeneous constraints on LAGRANGE >> (or HIERARCHIC, or pretty much anything but CLOUGH) elements. I don't >> know if anyone's ever actually tested using it before, though. > > I don't expect to be using anything other than LAGRANGE for a while... > > After adding DofConstraintRow's to the dof map, what generally needs to be > (re)initialized? Is it a good idea to call EquationSystems::reinit()? The difficulty is in how you add the rows  it needs to be done in user_constrain(), which gets called by EquationSystems::reinit(). If you do it any earlier your rows are probably going to be deleted when we construct the geometric/periodic constraints, and if you do it any later then your rows might not be taken into account properly during expansion of constraints. > I was certainly surprised... I'll try to verify that this is the only > difference between the working and nonworking cases and see if I can put > together a simple example. Thanks,  Roy 
From: Roy Stogner <roystgnr@ic...>  20100428 16:49:32

On Wed, 28 Apr 2010, liang cheng wrote: > I made a 3D finite deformation elasticity input code and got the executable > file, while the logical error displayed > in the screen when running it. > > Assertion `!_points.empty()' failed. > [0] /home/liang/Desktop/libmesh/include/quadrature/quadrature.h, line 116, This usually means that the quadrature rule wasn't properly reinitialized on an element. Did you attach it to your FE object, and did you reinit() the latter?  Roy 
From: liang cheng <goeasyon@gm...>  20100428 16:33:55

Hi developers and users, I made a 3D finite deformation elasticity input code and got the executable file, while the logical error displayed in the screen when running it. The Newmark solver is used for time steps and no Newton solver used because I would make the code simpler at beginning of trail. This local errors makes me wonder, I have no any clue to modify the code, does it lie in the main or assemble part, or even in the boundary condition? The assemble part is quit straightforward based on the Zienkiewicz's, the [K] martix follows the classic steps involving the displacement, deformation gradient, push forward the dphi[][], neoHookean algorithm and elastic tensor c_ijkl. The mass matrix [M] is set as lumped, which contains only diagonal components. No damping and body force here. The whole process looks fine but still get this error, does any one could give me some advices to solve it? Thanks in advance! Liang ************************** error ************************************************** liang@...:~/Desktop/libmesh/libmesh_liang/dld3d_v5$ ./dld3d_v5dbg Running ./dld3d_v5dbg Mesh Information: mesh_dimension()=3 spatial_dimension()=3 n_nodes()=2013 n_local_nodes()=2013 n_elem()=1200 n_local_elem()=1200 n_active_elem()=1200 n_subdomains()=1 n_processors()=1 processor_id()=0 EquationSystems n_systems()=1 System "Disp_finite_deform_3d" Type "Newmark"Error in libMesh internal logic Variables="u" "v" "w" Finite Element Types="LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00" Infinite Element Mapping="CARTESIAN" "CARTESIAN" "CARTESIAN" Approximation Orders="FIRST", "THIRD" "FIRST", "THIRD" "FIRST", "THIRD" n_dofs()=6039 n_local_dofs()=6039 n_constrained_dofs()=0 n_vectors()=9 Assertion `!_points.empty()' failed. [0] /home/liang/Desktop/libmesh/include/quadrature/quadrature.h, line 116, compiled Apr 28 2010 at 12:08:50 terminate called after throwing an instance of 'libMesh::LogicError' what(): Error in libMesh internal logic Aborted (core dumped) 
From: Jed Brown <jed@59...>  20100428 14:52:19

On Wed, 28 Apr 2010 12:06:47 +0200, Lorenzo Botti <bottilorenzo@...> wrote: > If my variables are u v w p (3D stokes equation) and I have an equal order > approximation > Can I just do > > pc_fieldsplit_0_fields 0,1,2 > pc_fieldsplit_1_fields 3 > pc_fieldsplit_type schur > pc_filedsplit_block_size <n_dofs_per_variable> Yes, and you can control the solver for the velocity block under fieldsplit_0_{ksp,pc} and for the Schur complement under fieldsplit_1_{ksp,pc}. Note that with some discretizations, like stabilized Q1Q1, you may get better results using an unsplit method (multigrid on the coupled system). This doesn't work with mixed methods unless you define custom restrictions and smoothing operators. Jed 
From: Lorenzo Botti <bottilorenzo@gm...>  20100428 10:06:54

Sorry to bother you, just to understand if I got it. If your variables are interlaced like > > u0,v0,w0,u1,v1,w1,... > > Then you can (a) set a block size for the matrix with MatSetBlockSize, > (b) use BAIJ in which case the block size is built in, or (c) use > pc_fieldsplit_block_size. Then to define the splits, use > > pc_fieldsplit_0_fields 0,3 pc_fieldsplit_1_fields 1 > pc_fieldsplit_2_fields 2,4,5 > > to split six fields into three splits of sizes 2,1,3. > > > If my variables are u v w p (3D stokes equation) and I have an equal order approximation Can I just do pc_fieldsplit_0_fields 0,1,2 pc_fieldsplit_1_fields 3 pc_fieldsplit_type schur pc_filedsplit_block_size <n_dofs_per_variable> Thank you very much Lorenzo 
From: Boyce Griffith <griffith@ci...>  20100428 03:29:25

On 4/27/10 10:20 PM, Roy Stogner wrote: > > On Tue, 27 Apr 2010, Boyce Griffith wrote: > >> I am trying to use the DofMap to impose Dirichlettype boundary >> conditions within an explicit nonlinear mechanics solver which is >> embedded in a larger fluidstructure interaction solver. I constrain >> the velocity and force DOFs at Dirichlet boundaries to be zero by adding >> a DofConstraintRow for each of the constrained DOFs. > > This should be a good method for homogeneous constraints on LAGRANGE > (or HIERARCHIC, or pretty much anything but CLOUGH) elements. I don't > know if anyone's ever actually tested using it before, though. I don't expect to be using anything other than LAGRANGE for a while... After adding DofConstraintRow's to the dof map, what generally needs to be (re)initialized? Is it a good idea to call EquationSystems::reinit()? >> This seems to work well as long as I do not also have subdomains in my >> mesh; however, as soon as I start labeling different subdomains, the dof >> constraints no longer impose the boundary conditions. Is this expected >> behavior? Is there something special that I should be doing in the case >> of multiple subdomains? > > This is very surprising  unless you're adding subdomainspecific > variables, in general the only part of the code that cares about > subdomains is the I/O (to read them in and write them out), the > refinement (to assign them to new child elements), and your own user > code! Could you verify this problem in a small example? I was certainly surprised... I'll try to verify that this is the only difference between the working and nonworking cases and see if I can put together a simple example. Thanks,  Boyce 
From: Roy Stogner <roystgnr@ic...>  20100428 02:21:02

On Tue, 27 Apr 2010, Boyce Griffith wrote: > I am trying to use the DofMap to impose Dirichlettype boundary > conditions within an explicit nonlinear mechanics solver which is > embedded in a larger fluidstructure interaction solver. I constrain > the velocity and force DOFs at Dirichlet boundaries to be zero by adding > a DofConstraintRow for each of the constrained DOFs. This should be a good method for homogeneous constraints on LAGRANGE (or HIERARCHIC, or pretty much anything but CLOUGH) elements. I don't know if anyone's ever actually tested using it before, though. > This seems to work well as long as I do not also have subdomains in my > mesh; however, as soon as I start labeling different subdomains, the dof > constraints no longer impose the boundary conditions. Is this expected > behavior? Is there something special that I should be doing in the case > of multiple subdomains? This is very surprising  unless you're adding subdomainspecific variables, in general the only part of the code that cares about subdomains is the I/O (to read them in and write them out), the refinement (to assign them to new child elements), and your own user code! Could you verify this problem in a small example? > By the way, subdomain_id seems to be the most direct way to assign > different material properties to different parts of the mesh. Are > there any alternatives which I should use instead? Subdomain ids are ideal unless you have a very large number of different materials or you have continuously varying material properties; in those cases you'd want to use a separate ExplicitSystem with properties as variable fields.  Roy 
From: Boyce Griffith <griffith@ci...>  20100428 00:21:53

Hi, Folks  I am trying to use the DofMap to impose Dirichlettype boundary conditions within an explicit nonlinear mechanics solver which is embedded in a larger fluidstructure interaction solver. I constrain the velocity and force DOFs at Dirichlet boundaries to be zero by adding a DofConstraintRow for each of the constrained DOFs. This seems to work well as long as I do not also have subdomains in my mesh; however, as soon as I start labeling different subdomains, the dof constraints no longer impose the boundary conditions. Is this expected behavior? Is there something special that I should be doing in the case of multiple subdomains? By the way, subdomain_id seems to be the most direct way to assign different material properties to different parts of the mesh. Are there any alternatives which I should use instead? Thanks!  Boyce 