From: Roy Stogner <roystgnr@ic...>  20110929 17:04:08

On Wed, 28 Sep 2011, Kirk, Benjamin (JSCEG311) wrote: > So am I right that the most pressing library change required is to > update the documentation? Sounds like. Personally *I* think the API is selfdocumenting: if the argument is nonconst then of course we're going to change it! ;) I guess the Doxygen comments really ought to specify why and to what it's being changed, though.  Roy 
From: David Andrs <David.A<ndrs@in...>  20110928 22:26:13
Attachments:
Message as HTML

/* $Id: ex4.C 2501 20071120 02:33:29Z benkirk $ */ /* 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 USA */ // <h1>Example 19  Solving the 2D Young Laplace Problem using nonlinear solvers</h1> // // This example shows how to use the NonlinearImplicitSystem class // to efficiently solve nonlinear problems in parallel. // // In nonlinear systems, we aim at finding x that satisfy R(x) = 0. // In nonlinear finite element analysis, the residual is typically // of the form R(x) = K(x)*x  f, with K(x) the system matrix and f // the "righthandside". The NonlinearImplicitSystem class expects // two callback functions to compute the residual R and its Jacobian // for the Newton iterations. Here, we just approximate // the true Jacobian by K(x). // // You can turn on preconditining of the matrix free system using the // jacobian by passing "pre" on the command line. Currently this only // work with Petsc so this isn't used by using "make run" // // This example also runs with the experimental Trilinos NOX solvers by specifying // the usetrilinos command line argument. // C++ include files that we need #include <iostream> #include <algorithm> #include <cmath> // Various include files needed for the mesh & solver functionality. #include "libmesh.h" #include "mesh.h" #include "mesh_refinement.h" #include "exodusII_io.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 "elem.h" #include "string_to_enum.h" #include "getpot.h" #include "dense_submatrix.h" #include "dense_subvector.h" #include "mesh_generation.h" #include "coupling_matrix.h" #include "boundary_info.h" #include "mesh_refinement.h" #include "error_vector.h" #include "kelly_error_estimator.h" // The nonlinear solver and system we will be using #include "nonlinear_solver.h" #include "nonlinear_implicit_system.h" // Necessary for programmatically setting petsc options #ifdef LIBMESH_HAVE_PETSC #include <petsc.h> #endif // Bring in everything from the libMesh namespace using namespace libMesh; // A reference to our equation system EquationSystems *_equation_system = NULL; // Lets define the physical parameters of the equation const Real kappa = 1.; const Real sigma = 0.2; // This function computes the Jacobian K(x) void compute_jacobian (const NumericVector<Number>& soln, SparseMatrix<Number>& jacobian, NonlinearImplicitSystem& sys) { // Get a reference to the equation system. EquationSystems &es = *_equation_system; // Get a constant reference to the mesh object. const MeshBase& mesh = es.get_mesh(); // The dimension that we are running const unsigned int dim = mesh.mesh_dimension(); // Get a reference to the NonlinearImplicitSystem we are solving NonlinearImplicitSystem& system = es.get_system<NonlinearImplicitSystem>("LaplaceYoung"); const unsigned int u_var = system.variable_number ("u"); const unsigned int v_var = system.variable_number ("v"); // A reference to the \p DofMap object for this system. The \p DofMap // object handles the index translation from node and element numbers // to degree of freedom numbers. We will talk more about the \p DofMap // in future examples. const DofMap& dof_map = system.get_dof_map(); // Get a constant reference to the Finite Element type // for the first (and only) variable in the system. FEType fe_type = dof_map.variable_type(0); // Build a Finite Element object of the specified type. Since the // \p FEBase::build() member dynamically creates memory we will // store the object as an \p AutoPtr<FEBase>. This can be thought // of as a pointer that will clean up after itself. AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); // A 5th order Gauss quadrature rule for numerical integration. QGauss qrule (dim, FIFTH); // Tell the finite element object to use our quadrature rule. fe>attach_quadrature_rule (&qrule); // Here we define some references to cellspecific data that // will be used to assemble the linear system. // We begin with the element Jacobian * quadrature weight at each // integration point. const std::vector<Real>& JxW = fe>get_JxW(); // The element shape functions evaluated at the quadrature points. const std::vector<std::vector<Real> >& phi = fe>get_phi(); // The element shape function gradients evaluated at the quadrature // points. const std::vector<std::vector<RealGradient> >& dphi = fe>get_dphi(); // Define data structures to contain the Jacobian element matrix. // Following basic finite element terminology we will denote these // "Ke". More detail is in example 3. // DenseMatrix<Number> Ke; // DenseMatrix<Number> Kuu(Ke), Kvv(Ke), Kuv(Ke); DenseMatrix<Number> Kuu, Kvv, Kuv; // This vector will hold the degree of freedom indices for // the element. These define where in the global system // the element degrees of freedom get mapped. std::vector<unsigned int> dof_indices; std::vector<unsigned int> dof_indices_u; std::vector<unsigned int> dof_indices_v; // Now we will loop over all the active elements in the mesh which // are local to this processor. // We will compute the element Jacobian contribution. MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); for ( ; el != end_el; ++el) { // Store a pointer to the element we are currently // working on. This allows for nicer syntax later. const Elem* elem = *el; // Get the degree of freedom indices for the // current element. These define where in the global // matrix and righthandside this element will // contribute to. 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 = dof_indices.size(); const unsigned int n_u_dofs = dof_indices_u.size(); const unsigned int n_v_dofs = dof_indices_v.size(); // Compute the elementspecific data for the current // element. This involves computing the location of the // quadrature points (q_point) and the shape functions // (phi, dphi) for the current element. fe>reinit (elem); // Zero the element Jacobian before // summing them. We use the resize member here because // the number of degrees of freedom might have changed from // the last element. Note that this will be the case if the // element type is different (i.e. the last element was a // triangle, now we are on a quadrilateral). // Ke.resize (dof_indices.size(), // dof_indices.size()); // Kuu.reposition( 0, 0, n_u_dofs, n_u_dofs); // Kuv.reposition(n_u_dofs, 0, n_u_dofs, n_v_dofs); // Kvv.reposition(n_u_dofs, n_u_dofs, n_v_dofs, n_v_dofs); Kuu.resize(n_u_dofs, n_u_dofs); Kuv.resize(n_u_dofs, n_v_dofs); Kvv.resize(n_v_dofs, n_v_dofs); // Now we will build the element Jacobian. This involves // a double loop to integrate the test funcions (i) against // the trial functions (j). Note that the Jacobian depends // on the current solution x, which we access using the soln // vector. // for (unsigned int qp=0; qp<qrule.n_points(); qp++) { Gradient grad_u; Gradient grad_v; for (unsigned int i=0; i<phi.size(); i++) { grad_u += dphi[i][qp]*soln(dof_indices_u[i]); grad_v += dphi[i][qp]*soln(dof_indices_v[i]); } for (unsigned int i=0; i<phi.size(); i++) for (unsigned int j=0; j<phi.size(); j++) { Kuu(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]); Kvv(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]); } } // dof_map.constrain_element_matrix (Ke, dof_indices); dof_map.constrain_element_matrix (Kuu, dof_indices_u, dof_indices_u); dof_map.constrain_element_matrix (Kuv, dof_indices_u, dof_indices_v); dof_map.constrain_element_matrix (Kvv, dof_indices_v, dof_indices_v); jacobian.add_matrix (Kuu, dof_indices_u, dof_indices_u); jacobian.add_matrix (Kuv, dof_indices_u, dof_indices_v); jacobian.add_matrix (Kvv, dof_indices_v, dof_indices_v); // Add the element matrix to the system Jacobian. // jacobian.add_matrix (Ke, dof_indices); } jacobian.close(); std::vector<int> rows; std::vector<unsigned int> nl; std::vector<short int> il; mesh.boundary_info>build_node_list(nl, il); for (int ii = 0; ii < nl.size(); ii++) { const Node & nd = mesh.node(nl[ii]); if (il[ii] == 1  il[ii] == 3) { unsigned int idx = nd.dof_number(system.number(), u_var, 0); rows.push_back(idx); } } jacobian.zero_rows(rows, 1); jacobian.close(); // That's it. } // Here we compute the residual R(x) = K(x)*x  f. The current solution // x is passed in the soln vector void compute_residual (const NumericVector<Number>& soln, NumericVector<Number>& residual, NonlinearImplicitSystem& sys) { EquationSystems &es = *_equation_system; // Get a constant reference to the mesh object. const MeshBase& mesh = es.get_mesh(); // The dimension that we are running const unsigned int dim = mesh.mesh_dimension(); libmesh_assert (dim == 2); // Get a reference to the NonlinearImplicitSystem we are solving NonlinearImplicitSystem& system = es.get_system<NonlinearImplicitSystem>("LaplaceYoung"); const unsigned int u_var = system.variable_number ("u"); const unsigned int v_var = system.variable_number ("v"); // A reference to the \p DofMap object for this system. The \p DofMap // object handles the index translation from node and element numbers // to degree of freedom numbers. We will talk more about the \p DofMap // in future examples. const DofMap& dof_map = system.get_dof_map(); // Get a constant reference to the Finite Element type // for the first (and only) variable in the system. FEType fe_type = dof_map.variable_type(0); // Build a Finite Element object of the specified type. Since the // \p FEBase::build() member dynamically creates memory we will // store the object as an \p AutoPtr<FEBase>. This can be thought // of as a pointer that will clean up after itself. AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); // A 5th order Gauss quadrature rule for numerical integration. QGauss qrule (dim, FIFTH); // Tell the finite element object to use our quadrature rule. fe>attach_quadrature_rule (&qrule); // Declare a special finite element object for // boundary integration. AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); // Boundary integration requires one quadraure rule, // with dimensionality one less than the dimensionality // of the element. QGauss qface(dim1, FIFTH); // Tell the finte element object to use our // quadrature rule. fe_face>attach_quadrature_rule (&qface); // Here we define some references to cellspecific data that // will be used to assemble the linear system. // We begin with the element Jacobian * quadrature weight at each // integration point. const std::vector<Real>& JxW = fe>get_JxW(); // The element shape functions evaluated at the quadrature points. const std::vector<std::vector<Real> >& phi = fe>get_phi(); // The element shape function gradients evaluated at the quadrature // points. const std::vector<std::vector<RealGradient> >& dphi = fe>get_dphi(); // Define data structures to contain the resdual contributions DenseVector<Number> Re; DenseSubVector<Number> Re_u(Re), Re_v(Re); // This vector will hold the degree of freedom indices for // the element. These define where in the global system // the element degrees of freedom get mapped. std::vector<unsigned int> dof_indices; std::vector<unsigned int> dof_indices_u, dof_indices_v; // Now we will loop over all the active elements in the mesh which // are local to this processor. // We will compute the element residual. residual.zero(); MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); for ( ; el != end_el; ++el) { // Store a pointer to the element we are currently // working on. This allows for nicer syntax later. const Elem* elem = *el; // Get the degree of freedom indices for the // current element. These define where in the global // matrix and righthandside this element will // contribute to. 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 = dof_indices.size(); const unsigned int n_u_dofs = dof_indices_u.size(); const unsigned int n_v_dofs = dof_indices_v.size(); // Compute the elementspecific data for the current // element. This involves computing the location of the // quadrature points (q_point) and the shape functions // (phi, dphi) for the current element. fe>reinit (elem); // We use the resize member here because // the number of degrees of freedom might have changed from // the last element. Note that this will be the case if the // element type is different (i.e. the last element was a // triangle, now we are on a quadrilateral). Re.resize (dof_indices.size()); Re_u.reposition(u_var * n_u_dofs, n_u_dofs); Re_v.reposition(v_var * n_u_dofs, n_v_dofs); // Now we will build the residual. This involves // the construction of the matrix K and multiplication of it // with the current solution x. We rearrange this into two loops: // In the first, we calculate only the contribution of // K_ij*x_j which is independent of the row i. In the second loops, // we multiply with the rowdependent part and add it to the element // residual. for (unsigned int qp=0; qp<qrule.n_points(); qp++) { Number u = 0; Gradient grad_u; Number v = 0; Gradient grad_v; for (unsigned int j=0; j<phi.size(); j++) { u += phi[j][qp]*soln(dof_indices_u[j]); grad_u += dphi[j][qp]*soln(dof_indices_u[j]); v += phi[j][qp]*soln(dof_indices_v[j]); grad_v += dphi[j][qp]*soln(dof_indices_v[j]); } for (unsigned int i=0; i<phi.size(); i++) { Re_u(i) += JxW[qp]*(dphi[i][qp]*grad_u); Re_v(i) += JxW[qp]*(dphi[i][qp]*grad_v); } } // At this point the interior element integration has // been completed. However, we have not yet addressed // boundary conditions. // // The following loops over the sides of the element. // // If the element has no neighbor on a side then that // // side MUST live on a boundary of the domain. // for (unsigned int side=0; side<elem>n_sides(); side++) // if (elem>neighbor(side) == NULL) // { // // The value of the shape functions at the quadrature // // points. // const std::vector<std::vector<Real> >& phi_face = fe_face>get_phi(); // // // The Jacobian * Quadrature Weight at the quadrature // // points on the face. // const std::vector<Real>& JxW_face = fe_face>get_JxW(); // // // Compute the shape function values on the element face. // fe_face>reinit(elem, side); // // // Loop over the face quadrature points for integration. // for (unsigned int qp=0; qp<qface.n_points(); qp++) // { // // This is the righthandside contribution (f), // // which has to be subtracted from the current residual // for (unsigned int i=0; i<phi_face.size(); i++) // Re_u(i) = JxW_face[qp]*sigma*phi_face[i][qp]; // } // } dof_map.constrain_element_vector (Re, dof_indices); residual.add_vector (Re, dof_indices); } residual.close(); std::vector<unsigned int> nl; std::vector<short int> il; mesh.boundary_info>build_node_list(nl, il); for (int ii = 0; ii < nl.size(); ii++) { const Node & nd = mesh.node(nl[ii]); if (il[ii] == 1) { unsigned int idx = nd.dof_number(system.number(), u_var, 0); residual.set(idx, soln(idx)  1.); } else if (il[ii] == 3) { unsigned int idx = nd.dof_number(system.number(), u_var, 0); residual.set(idx, soln(idx)  2.); } } residual.close(); // That's it. } // Begin the main program. int main (int argc, char** argv) { // Initialize libMesh and any dependent libaries, like in example 2. LibMeshInit init (argc, argv); #if !defined(LIBMESH_HAVE_PETSC) && !defined(LIBMESH_HAVE_TRILINOS) if (libMesh::processor_id() == 0) std::cerr << "ERROR: This example requires libMesh to be\n" << "compiled with nonlinear solver support from\n" << "PETSc or Trilinos!" << std::endl; return 0; #endif #ifndef LIBMESH_ENABLE_AMR if (libMesh::processor_id() == 0) std::cerr << "ERROR: This example requires libMesh to be\n" << "compiled with AMR support!" << std::endl; return 0; #else // Create a GetPot object to parse the command line GetPot command_line (argc, argv); // Skip this 2D example if libMesh was compiled as 1Donly. libmesh_example_assert(2 <= LIBMESH_DIM, "2D support"); // Create a mesh from file. Mesh mesh; MeshTools::Generation::build_square(mesh, 2, 2, 0, 1, 0, 1, QUAD4); mesh.boundary_info>build_node_list_from_side_list(); // Print information about the mesh to the screen. mesh.print_info(); // Create an equation systems object. EquationSystems equation_systems (mesh); _equation_system = &equation_systems; // Declare the system and its variables. // Creates a system named "LaplaceYoung" NonlinearImplicitSystem& system = equation_systems.add_system<NonlinearImplicitSystem> ("LaplaceYoung"); // Here we specify the tolerance for the nonlinear solver and // the maximum of nonlinear iterations. equation_systems.parameters.set<Real> ("nonlinear solver tolerance") = 1.e12; equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = 50; // Adds the variable "u" to "LaplaceYoung". "u" // will be approximated using secondorder approximation. system.add_variable("u", FIRST, LAGRANGE); system.add_variable("v", FIRST, LAGRANGE); CouplingMatrix cm(2); cm(0, 0) = 1; cm(1, 1) = 1; cm(0, 1) = 1; system.get_dof_map()._dof_coupling = &cm; // Give the system a pointer to the functions that update // the residual and Jacobian. system.nonlinear_solver>residual = compute_residual; system.nonlinear_solver>jacobian = compute_jacobian; // Initialize the data structures for the equation system. equation_systems.init(); // Prints information about the system to the screen. equation_systems.print_info(); MeshRefinement mesh_refinement (mesh); const unsigned int max_r_steps = 5; for (unsigned int r_step=0; r_step<max_r_steps; r_step++) { // Solve the system "LaplaceYoung", print the number of iterations // and final residual equation_systems.get_system("LaplaceYoung").solve(); // Print out final convergence information. This duplicates some // output from during the solve itself, but demonstrates another way // to get this information after the solve is complete. std::cout << "LaplaceYoung system solved at nonlinear iteration " << system.n_nonlinear_iterations() << " , final nonlinear residual norm: " << system.final_nonlinear_residual() << std::endl; if (r_step+1 != max_r_steps) { std::cout << " Refining the mesh..." << std::endl; ErrorVector error; KellyErrorEstimator error_estimator; error_estimator.estimate_error (system, error); mesh_refinement.refine_fraction() = 0.01; mesh_refinement.coarsen_fraction() = 0.0; mesh_refinement.max_h_level() = 5; mesh_refinement.flag_elements_by_error_fraction (error); mesh_refinement.refine_and_coarsen_elements(); equation_systems.reinit (); mesh.boundary_info>build_node_list_from_side_list(); } } #ifdef LIBMESH_HAVE_EXODUS_API // After solving the system write the solution ExodusII_IO (mesh).write_equation_systems ("out.e", equation_systems); #endif // #ifdef LIBMESH_HAVE_EXODUS_API #endif // #ifndef LIBMESH_ENABLE_AMR // All done. return 0; } 
From: Roy Stogner <roystgnr@ic...>  20110928 23:02:29

On Wed, 28 Sep 2011, David Andrs wrote: > We ran into an issue with our framework here at INL when forming full jacobian matrix as a preconditioner with AMR. What we hit > in our problem was the different size of the dense matrix (local jacobian contributions) and the dof indices array. I traced the > problem down to the level of DofMap::constrain_element_matrix. I found out that this method can change the dof indices array > (actually the method that changes the dof indices array is DofMap::build_constraint_matrix.) and that triggered an assert in > libMesh (attached is a modified ex19 that exhibits this behavior). This only happens when there are constrains involved. I quess > that no one bumped into this since people usually call this constrain method just once, but we call it multiple times and trying > to reuse those arrays. I've run into this before; my workaround was to copy the array before calling the constrain method. > So my question: Is this known behavior? And is this used somewhere else in libMesh, like we call that constrain_element_matrix > and we use this modified array for something later on? It's known behavior (in fact we've now got DofMap::constrain_nothing() whose sole purpose is this behavior), and IIRC I've talked to users with app code depending on this behavior, but I don't think it was originally *desired* behavior for its own sake; it's just that build_constraint_matrix() needs to expand that vector out recursively, so that method itself "reuses" the modified array. Whether that recursive use counts as "somewhere else" is a matter of perspective, but I think Ben made a fine design choice here: if build_constraint matrix did avoid modifying its input argument then it would have to keep both the old and new vector, and whenever the constraint recursion went N levels deep we'd have N different vectors floating around. Whereas even for our user codes where the original unexpanded vector does get needed again afterwards, we can get away with only one extra vector copy and one extra line of code to make the copy. On the other hand, I'm sure nobody's benchmarked both alternatives (or the alternatives where build_constraint_matrix uses std::set, or the alternative where a nonrecursive "userlevel" build_constraint_matrix method makes a vector copy and calls a recursive "librarylevel" method...), and I doubt it makes a noticeable performance difference. The gripping hand is that at this point I wouldn't want to change the old method in any case (see: "app code may depend on this behavior"). I don't think it's worth adding a new method either, unless there's some flaw I'm overlooking in my workaround?  Roy 
From: David Andrs <David.A<ndrs@in...>  20110929 14:49:59
Attachments:
Message as HTML

Roy Stogner <roystgnr@...> wrote on 09/28/2011 05:02:21 PM: > > On Wed, 28 Sep 2011, David Andrs wrote: > > > We ran into an issue with our framework here at INL when forming > full jacobian matrix as a preconditioner with AMR. What we hit > > in our problem was the different size of the dense matrix (local > jacobian contributions) and the dof indices array. I traced the > > problem down to the level of DofMap::constrain_element_matrix. I > found out that this method can change the dof indices array > > (actually the method that changes the dof indices array is > DofMap::build_constraint_matrix.) and that triggered an assert in > > libMesh (attached is a modified ex19 that exhibits this behavior). > This only happens when there are constrains involved. I quess > > that no one bumped into this since people usually call this > constrain method just once, but we call it multiple times and trying > > to reuse those arrays. > > I've run into this before; my workaround was to copy the array before > calling the constrain method. > > > So my question: Is this known behavior? And is this used somewhere > else in libMesh, like we call that constrain_element_matrix > > and we use this modified array for something later on? > > It's known behavior (in fact we've now got DofMap::constrain_nothing() > whose sole purpose is this behavior), and IIRC I've talked to users > with app code depending on this behavior, but I don't think it was > originally *desired* behavior for its own sake; it's just that > build_constraint_matrix() needs to expand that vector out recursively, > so that method itself "reuses" the modified array. > > Whether that recursive use counts as "somewhere else" is a matter of > perspective, but I think Ben made a fine design choice here: if > build_constraint matrix did avoid modifying its input argument then it > would have to keep both the old and new vector, and whenever the > constraint recursion went N levels deep we'd have N different vectors > floating around. > > Whereas even for our user codes where the original unexpanded vector > does get needed again afterwards, we can get away with only one extra > vector copy and one extra line of code to make the copy. > > On the other hand, I'm sure nobody's benchmarked both alternatives (or > the alternatives where build_constraint_matrix uses std::set, or the > alternative where a nonrecursive "userlevel" build_constraint_matrix > method makes a vector copy and calls a recursive "librarylevel" > method...), and I doubt it makes a noticeable performance difference. > > The gripping hand is that at this point I wouldn't want to change the > old method in any case (see: "app code may depend on this behavior"). > I don't think it's worth adding a new method either, unless there's > some flaw I'm overlooking in my workaround? >  > Roy Thanks, Roy, for the explanation. I wanted to ask before I take any action on the fix. The workaround you suggested above was one of the solutions we came up with.  David 
From: Kirk, Benjamin (JSCEG311) <benjamin.kirk1@na...>  20110928 23:51:44

(sorry to topreply, damn phone...) I also was originally planning for the most common case  in fact there are no constraints for the current element, hence no copy is required... So am I right that the most pressing library change required is to update the documentation? Ben On Sep 28, 2011, at 6:02 PM, "Roy Stogner"nc <roystgnr@...>bal wrote: > > On Wed, 28 Sep 2011, David Andrs wrote: > >> We ran into an issue with our framework here at INL when forming full jacobian matrix as a preconditioner with AMR. What we hit >> in our problem was the different size of the dense matrix (local jacobian contributions) and the dof indices array. I traced the >> problem down to the level of DofMap::constrain_element_matrix. I found out that this method can change the dof indices array >> (actually the method that changes the dof indices array is DofMap::build_constraint_matrix.) and that triggered an assert in >> libMesh (attached is a modified ex19 that exhibits this behavior). This only happens when there are constrains involved. I quess >> that no one bumped into this since people usually call this constrain method just once, but we call it multiple times and trying >> to reuse those arrays. > > I've run into this before; my workaround was to copy the array before > calling the constrain method. > >> So my question: Is this known behavior? And is this used somewhere else in libMesh, like we call that constrain_element_matrix >> and we use this modified array for something later on? > > It's known behavior (in fact we've now got DofMap::constrain_nothing() > whose sole purpose is this behavior), and IIRC I've talked to users > with app code depending on this behavior, but I don't think it was > originally *desired* behavior for its own sake; it's just that > build_constraint_matrix() needs to expand that vector out recursively, > so that method itself "reuses" the modified array. > > Whether that recursive use counts as "somewhere else" is a matter of > perspective, but I think Ben made a fine design choice here: if > build_constraint matrix did avoid modifying its input argument then it > would have to keep both the old and new vector, and whenever the > constraint recursion went N levels deep we'd have N different vectors > floating around. > > Whereas even for our user codes where the original unexpanded vector > does get needed again afterwards, we can get away with only one extra > vector copy and one extra line of code to make the copy. > > On the other hand, I'm sure nobody's benchmarked both alternatives (or > the alternatives where build_constraint_matrix uses std::set, or the > alternative where a nonrecursive "userlevel" build_constraint_matrix > method makes a vector copy and calls a recursive "librarylevel" > method...), and I doubt it makes a noticeable performance difference. > > The gripping hand is that at this point I wouldn't want to change the > old method in any case (see: "app code may depend on this behavior"). > I don't think it's worth adding a new method either, unless there's > some flaw I'm overlooking in my workaround? >  > Roy > >  > All the data continuously generated in your IT infrastructure contains a > definitive record of customers, application performance, security > threats, fraudulent activity and more. Splunk takes this data and makes > sense of it. Business sense. IT sense. Common sense. > http://p.sf.net/sfu/splunkd2dcopy1 > _______________________________________________ > Libmeshdevel mailing list > Libmeshdevel@... > https://lists.sourceforge.net/lists/listinfo/libmeshdevel 
From: Roy Stogner <roystgnr@ic...>  20110929 17:04:08

On Wed, 28 Sep 2011, Kirk, Benjamin (JSCEG311) wrote: > So am I right that the most pressing library change required is to > update the documentation? Sounds like. Personally *I* think the API is selfdocumenting: if the argument is nonconst then of course we're going to change it! ;) I guess the Doxygen comments really ought to specify why and to what it's being changed, though.  Roy 