Index: include/solvers/adaptive_time_solver.h =================================================================== --- include/solvers/adaptive_time_solver.h (revision 3205) +++ include/solvers/adaptive_time_solver.h (working copy) @@ -81,12 +81,14 @@ /** * This method is passed on to the core_time_solver */ - virtual bool element_residual (bool get_jacobian); + virtual bool element_residual (bool get_jacobian, + DiffContext&); /** * This method is passed on to the core_time_solver */ - virtual bool side_residual (bool get_jacobian); + virtual bool side_residual (bool get_jacobian, + DiffContext&); /** * An implicit linear or nonlinear solver to use at each timestep. Index: include/solvers/euler_solver.h =================================================================== --- include/solvers/euler_solver.h (revision 3205) +++ include/solvers/euler_solver.h (working copy) @@ -71,7 +71,8 @@ * to build a full residual on an element. What combination * it uses will depend on theta. */ - virtual bool element_residual (bool get_jacobian); + virtual bool element_residual (bool request_jacobian, + DiffContext&); /** * This method uses the DifferentiableSystem's @@ -79,7 +80,8 @@ * to build a full residual on an element's side. * What combination it uses will depend on theta. */ - virtual bool side_residual (bool get_jacobian); + virtual bool side_residual (bool request_jacobian, + DiffContext&); /** * The value for the theta method to employ: 1.0 corresponds Index: include/solvers/fem_system.h =================================================================== --- include/solvers/fem_system.h (revision 3205) +++ include/solvers/fem_system.h (working copy) @@ -27,17 +27,10 @@ // Local Includes #include "diff_system.h" -#include "elem.h" +#include "fem_context.h" -#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES -#include "tensor_value.h" -#endif - // Forward Declarations -class FEBase; -class QBase; -template class NumericVector; /** * This class provides a specific system class. It aims @@ -173,7 +166,8 @@ * formulations will require reimplementing eulerian_residual() * entirely. */ - virtual bool eulerian_residual (bool request_jacobian); + virtual bool eulerian_residual (bool request_jacobian, + DiffContext &context); /** * Adds a mass vector contribution on \p elem to elem_residual. @@ -188,9 +182,28 @@ * mass matrix (e.g. for divergence-free elements or mass lumping) * requires reimplementing mass_residual(). */ - virtual bool mass_residual (bool request_jacobian); + virtual bool mass_residual (bool request_jacobian, + DiffContext &context); /** + * Builds a FEMContext object with enough information to do + * evaluations on each element. + * + * For most problems, the default FEMSystem implementation is correct; users + * who subclass FEMContext will need to also reimplement this method to build + * it. + */ + virtual AutoPtr build_context(); + + /* + * Prepares the result of a build_context() call for use. + * + * Most FEMSystem-based problems will need to reimplement this in order to + * call FE::get_*() as their particular physics requires. + */ + virtual void init_context(DiffContext &); + + /** * Runs a postprocessing loop over all elements, and if * \p postprocess_sides is true over all sides. */ @@ -216,102 +229,6 @@ int extra_quadrature_order; /** - * Returns the value of the solution variable \p var at the quadrature - * point \p qp on the current element interior - */ - Number interior_value(unsigned int var, unsigned int qp); - - /** - * Returns the value of the solution variable \p var at the quadrature - * point \p qp on the current element side - */ - Number side_value(unsigned int var, unsigned int qp); - - /** - * Returns the value of the solution variable \p var at the physical - * point \p p on the current element - */ - Number point_value(unsigned int var, Point &p); - - /** - * Returns the gradient of the solution variable \p var at the quadrature - * point \p qp on the current element interior - */ - Gradient interior_gradient(unsigned int var, unsigned int qp); - - /** - * Returns the gradient of the solution variable \p var at the quadrature - * point \p qp on the current element side - */ - Gradient side_gradient(unsigned int var, unsigned int qp); - -#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES - /** - * Returns the hessian of the solution variable \p var at the quadrature - * point \p qp on the current element interior - */ - Tensor interior_hessian(unsigned int var, unsigned int qp); - - /** - * Returns the hessian of the solution variable \p var at the quadrature - * point \p qp on the current element side - */ - Tensor side_hessian(unsigned int var, unsigned int qp); - -#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES - - /** - * Returns the value of the fixed_solution variable \p var at the quadrature - * point \p qp on the current element interior - */ - Number fixed_interior_value(unsigned int var, unsigned int qp); - - /** - * Returns the value of the fixed_solution variable \p var at the quadrature - * point \p qp on the current element side - */ - Number fixed_side_value(unsigned int var, unsigned int qp); - - /** - * Returns the value of the fixed_solution variable \p var at the physical - * point \p p on the current element - */ - Number fixed_point_value(unsigned int var, Point &p); - - /** - * Returns the gradient of the fixed_solution variable \p var at the quadrature - * point \p qp on the current element interior - */ - Gradient fixed_interior_gradient(unsigned int var, unsigned int qp); - - /** - * Returns the gradient of the fixed_solution variable \p var at the quadrature - * point \p qp on the current element side - */ - Gradient fixed_side_gradient(unsigned int var, unsigned int qp); - -#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES - /** - * Returns the hessian of the fixed_solution variable \p var at the quadrature - * point \p qp on the current element interior - */ - Tensor fixed_interior_hessian(unsigned int var, unsigned int qp); - - /** - * Returns the hessian of the fixed_solution variable \p var at the quadrature - * point \p qp on the current element side - */ - Tensor fixed_side_hessian(unsigned int var, unsigned int qp); - -#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES - - /** - * @returns \p "General". Helps in identifying - * the system type in an equation system file. - */ -// virtual std::string system_type () const { return "PDE"; } - - /** * If calculating numeric jacobians is required, the FEMSystem * will perturb each solution vector entry by numerical_jacobian_h * when calculating finite differences. @@ -333,22 +250,7 @@ */ Real verify_analytic_jacobians; - /** - * Reinitialize Elem and FE objects if necessary for integration at a new - * point in time: specifically, handle moving elements in moving mesh - * schemes. - */ - virtual void elem_reinit(Real theta); - - /** - * Reinitialize Elem and side FE objects if necessary for integration at a - * new point in time: specifically, handle moving elements in moving mesh - * schemes. - */ - virtual void elem_side_reinit(Real theta); - protected: - /** * Initializes the member data fields associated with * the system, so that, e.g., \p assemble() may be used. @@ -356,91 +258,32 @@ virtual void init_data (); /** - * Clear data pointers associated with this system. - */ - void clear_fem_ptrs (); - - /** - * Uses the coordinate data specified by mesh_*_position configuration - * to set the geometry of \p elem to the value it would take after a fraction - * \p theta of a timestep. - */ - void elem_position_set(Real theta); - - /** - * Uses the geometry of \p elem to set the coordinate data specified - * by mesh_*_position configuration. - */ - void elem_position_get(); - - /** - * Reinitializes interior FE objects on the current geometric element - */ - void elem_fe_reinit(); - - /** - * Reinitializes side FE objects on the current geometric element - */ - void elem_side_fe_reinit(); - - - /** * Syntax sugar to make numerical_jacobian() declaration easier. */ - typedef bool (TimeSolver::*TimeSolverResPtr)(bool); + typedef bool (TimeSolver::*TimeSolverResPtr)(bool, DiffContext&); /** * Uses the results of multiple \p res calls * to numerically differentiate the corresponding jacobian. */ - void numerical_jacobian (TimeSolverResPtr res); + void numerical_jacobian (TimeSolverResPtr res, + FEMContext &context); /** * Uses the results of multiple element_residual() calls * to numerically differentiate the corresponding jacobian * on an element. */ - void numerical_elem_jacobian (); + void numerical_elem_jacobian (FEMContext &context); /** * Uses the results of multiple side_residual() calls * to numerically differentiate the corresponding jacobian * on an element's side. */ - void numerical_side_jacobian (); + void numerical_side_jacobian (FEMContext &context); /** - * Finite element objects for each variable's interior and sides. - */ - std::map element_fe; - std::map side_fe; - - /** - * Pointers to the same finite element objects, but indexed - * by variable number - */ - std::vector element_fe_var; - std::vector side_fe_var; - - /** - * Quadrature rules for element interior and sides. - * The FEM system will try to find a quadrature rule that - * correctly integrates all variables - */ - QBase *element_qrule; - QBase *side_qrule; - - /** - * Current element for element_* to examine - */ - Elem *elem; - - /** - * Current element for side_* to examine - */ - unsigned int side; - - /** * System and variables from which to acquire moving mesh information */ unsigned int _mesh_sys, _mesh_x_var, _mesh_y_var, _mesh_z_var; Index: include/solvers/fem_context.h =================================================================== --- include/solvers/fem_context.h (revision 0) +++ include/solvers/fem_context.h (revision 0) @@ -0,0 +1,250 @@ + +// $Id: fem_system.h 3103 2008-10-16 21:09:35Z roystgnr $ + +// The libMesh Finite Element Library. +// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// 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 02111-1307 USA + + + +#ifndef __fem_context_h__ +#define __fem_context_h__ + +// C++ includes + +// Local Includes +#include "diff_context.h" +#include "elem.h" + +#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES +#include "tensor_value.h" +#endif + +// Forward Declarations + +class FEBase; +class FEMSystem; +class QBase; +template class NumericVector; + +/** + * This class provides all data required for a physics package + * (e.g. an FEMSystem subclass) to perform local element residual + * and jacobian integrations. + * + * This class is part of the new DifferentiableSystem framework, + * which is still experimental. Users of this framework should + * beware of bugs and future API changes. + * + * @author Roy H. Stogner 2009 + */ + +// ------------------------------------------------------------ +// FEMContext class definition + +class FEMContext : public DiffContext +{ +public: + + /** + * Constructor. Allocates some but fills no data structures. + */ + FEMContext (const FEMSystem &sys); + + /** + * Destructor. + */ + virtual ~FEMContext (); + + /** + * Returns the value of the solution variable \p var at the quadrature + * point \p qp on the current element interior + */ + Number interior_value(unsigned int var, unsigned int qp); + + /** + * Returns the value of the solution variable \p var at the quadrature + * point \p qp on the current element side + */ + Number side_value(unsigned int var, unsigned int qp); + + /** + * Returns the value of the solution variable \p var at the physical + * point \p p on the current element + */ + Number point_value(unsigned int var, Point &p); + + /** + * Returns the gradient of the solution variable \p var at the quadrature + * point \p qp on the current element interior + */ + Gradient interior_gradient(unsigned int var, unsigned int qp); + + /** + * Returns the gradient of the solution variable \p var at the quadrature + * point \p qp on the current element side + */ + Gradient side_gradient(unsigned int var, unsigned int qp); + +#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES + /** + * Returns the hessian of the solution variable \p var at the quadrature + * point \p qp on the current element interior + */ + Tensor interior_hessian(unsigned int var, unsigned int qp); + + /** + * Returns the hessian of the solution variable \p var at the quadrature + * point \p qp on the current element side + */ + Tensor side_hessian(unsigned int var, unsigned int qp); + +#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES + + /** + * Returns the value of the fixed_solution variable \p var at the quadrature + * point \p qp on the current element interior + */ + Number fixed_interior_value(unsigned int var, unsigned int qp); + + /** + * Returns the value of the fixed_solution variable \p var at the quadrature + * point \p qp on the current element side + */ + Number fixed_side_value(unsigned int var, unsigned int qp); + + /** + * Returns the value of the fixed_solution variable \p var at the physical + * point \p p on the current element + */ + Number fixed_point_value(unsigned int var, Point &p); + + /** + * Returns the gradient of the fixed_solution variable \p var at the quadrature + * point \p qp on the current element interior + */ + Gradient fixed_interior_gradient(unsigned int var, unsigned int qp); + + /** + * Returns the gradient of the fixed_solution variable \p var at the quadrature + * point \p qp on the current element side + */ + Gradient fixed_side_gradient(unsigned int var, unsigned int qp); + +#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES + /** + * Returns the hessian of the fixed_solution variable \p var at the quadrature + * point \p qp on the current element interior + */ + Tensor fixed_interior_hessian(unsigned int var, unsigned int qp); + + /** + * Returns the hessian of the fixed_solution variable \p var at the quadrature + * point \p qp on the current element side + */ + Tensor fixed_side_hessian(unsigned int var, unsigned int qp); + +#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES + + /** + * Reinitialize all my FEM context data on a given + * element for the given system + */ + virtual void reinit(const FEMSystem&, Elem*); + +// should be protected: + /** + * Reinitialize Elem and FE objects if necessary for integration at a new + * point in time: specifically, handle moving elements in moving mesh + * schemes. + */ + virtual void elem_reinit(Real theta); + + /** + * Reinitialize Elem and side FE objects if necessary for integration at a + * new point in time: specifically, handle moving elements in moving mesh + * schemes. + */ + virtual void elem_side_reinit(Real theta); + + /** + * Reinitializes interior FE objects on the current geometric element + */ + void elem_fe_reinit(); + + /** + * Reinitializes side FE objects on the current geometric element + */ + void elem_side_fe_reinit(); + +// should be protected?: + /** + * Finite element objects for each variable's interior and sides. + */ + std::map element_fe; + std::map side_fe; + + /** + * Pointers to the same finite element objects, but indexed + * by variable number + */ + std::vector element_fe_var; + std::vector side_fe_var; + + /** + * Quadrature rules for element interior and sides. + * The FEM system will try to find a quadrature rule that + * correctly integrates all variables + */ + QBase *element_qrule; + QBase *side_qrule; + + /** + * Uses the coordinate data specified by mesh_*_position configuration + * to set the geometry of \p elem to the value it would take after a fraction + * \p theta of a timestep. + */ + void elem_position_set(Real theta); + + /** + * Uses the geometry of \p elem to set the coordinate data specified + * by mesh_*_position configuration. + */ + void elem_position_get(); + + /** + * System and variables from which to acquire moving mesh information + */ + unsigned int _mesh_sys, _mesh_x_var, _mesh_y_var, _mesh_z_var; + + /** + * Current element for element_* to examine + */ + Elem *elem; + + /** + * Current side for side_* to examine + */ + unsigned char side; + + /** + * Cached dimension of elements in this mesh + */ + unsigned char dim; +}; + + +#endif Index: include/solvers/diff_context.h =================================================================== --- include/solvers/diff_context.h (revision 0) +++ include/solvers/diff_context.h (revision 0) @@ -0,0 +1,153 @@ + +// $Id: diff_system.h 3086 2008-10-07 20:28:08Z roystgnr $ + +// The libMesh Finite Element Library. +// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// 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 02111-1307 USA + + + +#ifndef __diff_context_h__ +#define __diff_context_h__ + +// C++ includes + +// Local Includes +#include "dense_matrix.h" +#include "dense_submatrix.h" +#include "dense_subvector.h" +#include "dense_vector.h" + +// Forward declarations +class DifferentiableSystem; + +/** + * This class provides all data required for a physics package + * (e.g. a DifferentiableSystem subclass) to perform local element + * residual and jacobian integrations. + * + * This class is part of the new DifferentiableSystem framework, + * which is still experimental. Users of this framework should + * beware of bugs and future API changes. + * + * @author Roy H. Stogner 2009 + */ + +// ------------------------------------------------------------ +// DifferentiableSystem class definition + +class DiffContext +{ +public: + + /** + * Constructor. Optionally initializes required + * data structures. + */ + DiffContext (const DifferentiableSystem &); + + /** + * Destructor. + */ + virtual ~DiffContext (); + + /** + * If \p postprocess_sides is true (it is false by default), the + * postprocessing loop will loop over all sides as well as all elements. + */ + bool postprocess_sides; + + /** + * Gives derived classes the opportunity to reinitialize data (FE objects in + * FEMSystem, for example) needed for an interior integration at a new point + * within a timestep + */ + virtual void elem_reinit(Real) {} + + /** + * Gives derived classes the opportunity to reinitialize data needed for a + * side integration at a new point within a timestep + */ + virtual void elem_side_reinit(Real) {} + + /** + * For time-dependent problems, this is the time t for which the current + * nonlinear_solution is defined. + * FIXME - this needs to be tweaked mid-timestep by all transient solvers! + */ + Real time; + + /** + * Element by element components of nonlinear_solution + * as adjusted by a time_solver + */ + DenseVector elem_solution; + std::vector *> elem_subsolutions; + + /** + * Element by element components of nonlinear_solution + * at a fixed point in a timestep, for optional use by e.g. + * stabilized methods + */ + DenseVector elem_fixed_solution; + std::vector *> elem_fixed_subsolutions; + + /** + * A boolean to be set to true by the library whenever a jacobian + * evaluation is being requested. + */ + bool request_jacobian; + + /** + * The derivative of elem_solution with respect to the nonlinear solution, + * for use by systems constructing jacobians with elem_fixed_solution + * based methods + */ + Real elem_solution_derivative; + + /** + * The derivative of elem_fixed_solution with respect to the nonlinear + * solution, for use by systems constructing jacobians with + * elem_fixed_solution based methods + */ + Real fixed_solution_derivative; + + /** + * Element residual vector + */ + DenseVector elem_residual; + + /** + * Element jacobian: derivatives of elem_residual with respect to + * elem_solution + */ + DenseMatrix elem_jacobian; + + /** + * Element residual subvectors and Jacobian submatrices + */ + std::vector *> elem_subresiduals; + std::vector *> > elem_subjacobians; + + /** + * Global Degree of freedom index lists + */ + std::vector dof_indices; + std::vector > dof_indices_var; +}; + + +#endif Index: include/solvers/steady_solver.h =================================================================== --- include/solvers/steady_solver.h (revision 3205) +++ include/solvers/steady_solver.h (working copy) @@ -29,8 +29,9 @@ #include "time_solver.h" // Forward Declarations +class DiffContext; +class DifferentiableSystem; class TimeSolver; -class DifferentiableSystem; /** * This class implements a TimeSolver which does a single @@ -80,14 +81,16 @@ * element_time_derivative() and element_constraint() * to build a full residual/jacobian on an element. */ - virtual bool element_residual (bool get_jacobian); + virtual bool element_residual (bool request_jacobian, + DiffContext &); /** * This method uses the DifferentiableSystem's * side_time_derivative() and side_constraint() * to build a full residual/jacobian on an element's side. */ - virtual bool side_residual (bool get_jacobian); + virtual bool side_residual (bool request_jacobian, + DiffContext &); /** * Nominally computes the size of the difference between Index: include/solvers/euler2_solver.h =================================================================== --- include/solvers/euler2_solver.h (revision 3205) +++ include/solvers/euler2_solver.h (working copy) @@ -76,7 +76,8 @@ * to build a full residual on an element. What combination * it uses will depend on theta. */ - virtual bool element_residual (bool get_jacobian); + virtual bool element_residual (bool request_jacobian, + DiffContext&); /** * This method uses the DifferentiableSystem's @@ -84,7 +85,8 @@ * to build a full residual on an element's side. * What combination it uses will depend on theta. */ - virtual bool side_residual (bool get_jacobian); + virtual bool side_residual (bool request_jacobian, + DiffContext&); /** * The value for the theta method to employ: 1.0 corresponds Index: include/solvers/diff_system.h =================================================================== --- include/solvers/diff_system.h (revision 3205) +++ include/solvers/diff_system.h (working copy) @@ -30,6 +30,7 @@ #include "dense_submatrix.h" #include "dense_subvector.h" #include "dense_vector.h" +#include "diff_context.h" #include "linear_implicit_system.h" #include "transient_system.h" @@ -115,7 +116,8 @@ * examine u = elem_solution and add (F(u), phi_i) to * elem_residual. */ - virtual bool element_time_derivative (bool request_jacobian) { + virtual bool element_time_derivative (bool request_jacobian, + DiffContext &) { return request_jacobian; } @@ -131,7 +133,8 @@ * examine u = elem_solution and add (G(u), phi_i) to * elem_residual. */ - virtual bool element_constraint (bool request_jacobian) { + virtual bool element_constraint (bool request_jacobian, + DiffContext &) { return request_jacobian; } @@ -154,7 +157,8 @@ * Users may need to reimplement this for their particular PDE * depending on the boundary conditions. */ - virtual bool side_time_derivative (bool request_jacobian) { + virtual bool side_time_derivative (bool request_jacobian, + DiffContext &) { return request_jacobian; } @@ -169,7 +173,8 @@ * Users may need to reimplement this for their particular PDE * depending on the boundary conditions. */ - virtual bool side_constraint (bool request_jacobian) { + virtual bool side_constraint (bool request_jacobian, + DiffContext &) { return request_jacobian; } @@ -188,13 +193,13 @@ /** * Does any work that needs to be done on \p elem in a postprocessing loop. */ - virtual void element_postprocess () {} + virtual void element_postprocess (DiffContext &) {} /** * Does any work that needs to be done on \p side of \p elem in a * postprocessing loop. */ - virtual void side_postprocess () {} + virtual void side_postprocess (DiffContext &) {} /** * Tells the DiffSystem that variable var is evolving with @@ -220,7 +225,8 @@ * The library provides a basic implementation in * FEMSystem::eulerian_residual() */ - virtual bool eulerian_residual (bool request_jacobian) { + virtual bool eulerian_residual (bool request_jacobian, + DiffContext &) { return request_jacobian; } @@ -235,7 +241,8 @@ * FEMSystem::mass_residual; few users will need to reimplement * this themselves. */ - virtual bool mass_residual (bool request_jacobian) { + virtual bool mass_residual (bool request_jacobian, + DiffContext &) { return request_jacobian; } @@ -251,23 +258,29 @@ * is correct; users with boundary conditions including time * derivatives may need to reimplement this themselves. */ - virtual bool side_mass_residual (bool request_jacobian) { + virtual bool side_mass_residual (bool request_jacobian, + DiffContext &) { return request_jacobian; } /** - * Gives derived classes the opportunity to reinitialize data (FE objects in - * FEMSystem, for example) needed for an interior integration at a new point - * within a timestep + * Builds a DiffContext object with enough information to do + * evaluations on each element. + * + * For most problems, the default "Let FEMSystem build an FEMContext + * reimplementation is correct; users who subclass FEMContext will need to + * also reimplement this method to build it. */ - virtual void elem_reinit(Real) {} + virtual AutoPtr build_context(); - /** - * Gives derived classes the opportunity to reinitialize data needed for a - * side integration at a new point within a timestep + /* + * Prepares the result of a build_context() call for use. + * + * Most FEMSystem-based problems will need to reimplement this in order to + * call FE::get_*() as their particular physics requires. */ - virtual void elem_side_reinit(Real) {} - + virtual void init_context(DiffContext &) {} + /** * Invokes the solver associated with the system. For steady state * solvers, this will find a root x where F(x) = 0. For transient @@ -276,20 +289,16 @@ virtual void solve (); /** - * @returns \p "Differentiable". Helps in identifying - * the system type in an equation system file. - */ -// virtual std::string system_type () const { return "Differentiable"; } - - /** * A pointer to the solver object we're going to use. * This must be instantiated by the user before solving! */ AutoPtr time_solver; /** - * For time-dependent problems, this is the time t for which the current - * nonlinear_solution is defined. + * For time-dependent problems, this is the time t at the beginning of + * the current timestep. But do *not* access this time during an + * assembly! Use the DiffContext::time value instead to get correct + * results. */ Real time; @@ -337,21 +346,6 @@ bool print_element_jacobians; /** - * Element by element components of nonlinear_solution - * as adjusted by a time_solver - */ - DenseVector elem_solution; - std::vector *> elem_subsolutions; - - /** - * Element by element components of nonlinear_solution - * at a fixed point in a timestep, for optional use by e.g. - * stabilized methods - */ - DenseVector elem_fixed_solution; - std::vector *> elem_fixed_subsolutions; - - /** * A boolean to be set to true by systems using elem_fixed_solution, * so that the library code will create it. * False by default. @@ -360,43 +354,6 @@ */ bool use_fixed_solution; - /** - * The derivative of elem_solution with respect to the nonlinear solution, - * for use by systems constructing jacobians with elem_fixed_solution - * based methods - */ - Real elem_solution_derivative; - - /** - * The derivative of elem_fixed_solution with respect to the nonlinear - * solution, for use by systems constructing jacobians with - * elem_fixed_solution based methods - */ - Real fixed_solution_derivative; - - /** - * Element residual vector - */ - DenseVector elem_residual; - - /** - * Element jacobian: derivatives of elem_residual with respect to - * elem_solution - */ - DenseMatrix elem_jacobian; - - /** - * Element residual subvectors and Jacobian submatrices - */ - std::vector *> elem_subresiduals; - std::vector *> > elem_subjacobians; - - /** - * Global Degree of freedom index lists - */ - std::vector dof_indices; - std::vector > dof_indices_var; - protected: /** * Initializes the member data fields associated with @@ -409,11 +366,6 @@ * in time and which are just constraints */ std::vector _time_evolving; - - /** - * Clear data pointers associated with this system. - */ - void clear_diff_ptrs (); }; Index: include/solvers/time_solver.h =================================================================== --- include/solvers/time_solver.h (revision 3205) +++ include/solvers/time_solver.h (working copy) @@ -32,6 +32,7 @@ #include "reference_counted_object.h" // Forward Declarations +class DiffContext; class DiffSolver; class TimeSolver; class DifferentiableSystem; @@ -106,7 +107,8 @@ * it uses will depend on the type of solver. See * the subclasses for more details. */ - virtual bool element_residual (bool get_jacobian) = 0; + virtual bool element_residual (bool request_jacobian, + DiffContext &) = 0; /** * This method uses the DifferentiableSystem's @@ -115,7 +117,8 @@ * What combination it uses will depend on the type * of solver. See the subclasses for more details. */ - virtual bool side_residual (bool get_jacobian) = 0; + virtual bool side_residual (bool request_jacobian, + DiffContext &) = 0; /** * This method is for subclasses or users to override Index: src/solvers/diff_system.C =================================================================== --- src/solvers/diff_system.C (revision 3205) +++ src/solvers/diff_system.C (working copy) @@ -42,9 +42,7 @@ print_jacobian_norms(false), print_jacobians(false), print_element_jacobians(false), - use_fixed_solution(false), - elem_solution_derivative(1.), - fixed_solution_derivative(0.) + use_fixed_solution(false) { untested(); } @@ -60,8 +58,6 @@ void DifferentiableSystem::clear () { - this->clear_diff_ptrs(); - _time_evolving.resize(0); use_fixed_solution = false; @@ -69,31 +65,6 @@ -void DifferentiableSystem::clear_diff_ptrs () -{ - for (unsigned int i=0; i != elem_subsolutions.size(); ++i) - { - delete elem_subsolutions[i]; - delete elem_subresiduals[i]; - - if (use_fixed_solution) - { - delete elem_fixed_subsolutions[i]; - } - - for (unsigned int j=0; j != elem_subjacobians[i].size(); ++j) - delete elem_subjacobians[i][j]; - } - elem_subsolutions.resize(0); - elem_subresiduals.resize(0); - elem_subjacobians.resize(0); - - if (use_fixed_solution) - elem_fixed_subsolutions.resize(0); -} - - - void DifferentiableSystem::reinit () { Parent::reinit(); @@ -105,19 +76,7 @@ void DifferentiableSystem::init_data () { - // First, allocate a vector for the iterate in our quasi_Newton - // solver - - // Update - we now just use the System - // solution vector if the solver doesn't need an extra vector - - // We don't want to project more solutions than we have to -// this->add_vector("_nonlinear_solution", false); -// this->add_vector("_nonlinear_solution"); -// this->project_solution_on_reinit() = false; - - // Next, give us flags for every variable that might be time - // evolving + // give us flags for every variable that might be time evolving _time_evolving.resize(this->n_vars(), false); // Do any initialization our solvers need @@ -126,43 +85,13 @@ // Next initialize ImplicitSystem data Parent::init_data(); +} - // Finally initialize solution/residual/jacobian data structures - unsigned int n_vars = this->n_vars(); - dof_indices_var.resize(n_vars); - // We may have already been initialized - this->clear_diff_ptrs(); - - elem_subsolutions.reserve(n_vars); - elem_subresiduals.reserve(n_vars); - elem_subjacobians.resize(n_vars); - - if (use_fixed_solution) - elem_fixed_subsolutions.reserve(n_vars); - - for (unsigned int i=0; i != n_vars; ++i) - { - elem_subsolutions.push_back(new DenseSubVector(elem_solution)); - elem_subresiduals.push_back(new DenseSubVector(elem_residual)); - elem_subjacobians[i].clear(); - elem_subjacobians[i].reserve(n_vars); - - if (use_fixed_solution) - { - elem_fixed_subsolutions.push_back - (new DenseSubVector(elem_fixed_solution)); - } - for (unsigned int j=0; j != n_vars; ++j) - { - elem_subjacobians[i].push_back - (new DenseSubMatrix(elem_jacobian)); - } - } - - // Resize the serial nonlinear solution for the current mesh - // current_local_nonlinear_solution->init (this->n_dofs()); +AutoPtr DifferentiableSystem::build_context () +{ + return AutoPtr(new DiffContext(*this)); } Index: src/solvers/adaptive_time_solver.C =================================================================== --- src/solvers/adaptive_time_solver.C (revision 3205) +++ src/solvers/adaptive_time_solver.C (working copy) @@ -108,20 +108,22 @@ -bool AdaptiveTimeSolver::element_residual (bool request_jacobian) +bool AdaptiveTimeSolver::element_residual (bool request_jacobian, + DiffContext &context) { libmesh_assert(core_time_solver.get()); - return core_time_solver->element_residual(request_jacobian); + return core_time_solver->element_residual(request_jacobian, context); } -bool AdaptiveTimeSolver::side_residual (bool request_jacobian) +bool AdaptiveTimeSolver::side_residual (bool request_jacobian, + DiffContext &context) { libmesh_assert(core_time_solver.get()); - return core_time_solver->side_residual(request_jacobian); + return core_time_solver->side_residual(request_jacobian, context); } Index: src/solvers/euler_solver.C =================================================================== --- src/solvers/euler_solver.C (revision 3205) +++ src/solvers/euler_solver.C (working copy) @@ -47,68 +47,69 @@ -bool EulerSolver::element_residual (bool request_jacobian) +bool EulerSolver::element_residual (bool request_jacobian, + DiffContext &context) { - unsigned int n_dofs = _system.elem_solution.size(); + unsigned int n_dofs = context.elem_solution.size(); // Local nonlinear solution at old timestep DenseVector old_elem_solution(n_dofs); for (unsigned int i=0; i != n_dofs; ++i) old_elem_solution(i) = - old_nonlinear_solution(_system.dof_indices[i]); + old_nonlinear_solution(context.dof_indices[i]); // Local nonlinear solution at time t_theta - DenseVector theta_solution(_system.elem_solution); + DenseVector theta_solution(context.elem_solution); theta_solution *= theta; theta_solution.add(1. - theta, old_elem_solution); // Technically the elem_solution_derivative is either theta // or -1.0 in this implementation, but we scale the former part // ourselves - _system.elem_solution_derivative = 1.0; + context.elem_solution_derivative = 1.0; // Technically the fixed_solution_derivative is always theta, // but we're scaling a whole jacobian by theta after these first // evaluations - _system.fixed_solution_derivative = 1.0; + context.fixed_solution_derivative = 1.0; // If a fixed solution is requested, we'll use theta_solution if (_system.use_fixed_solution) - _system.elem_fixed_solution = theta_solution; + context.elem_fixed_solution = theta_solution; // Move theta_->elem_, elem_->theta_ - _system.elem_solution.swap(theta_solution); + context.elem_solution.swap(theta_solution); // Move the mesh into place first if necessary - _system.elem_reinit(theta); + context.elem_reinit(theta); // We're going to compute just the change in elem_residual // (and possibly elem_jacobian), then add back the old values - DenseVector old_elem_residual(_system.elem_residual); + DenseVector old_elem_residual(context.elem_residual); DenseMatrix old_elem_jacobian; if (request_jacobian) { - old_elem_jacobian = _system.elem_jacobian; - _system.elem_jacobian.zero(); + old_elem_jacobian = context.elem_jacobian; + context.elem_jacobian.zero(); } - _system.elem_residual.zero(); + context.elem_residual.zero(); // Get the time derivative at t_theta bool jacobian_computed = - _system.element_time_derivative(request_jacobian); + _system.element_time_derivative(request_jacobian, context); // For a moving mesh problem we may need the pseudoconvection term too jacobian_computed = - _system.eulerian_residual(jacobian_computed) && jacobian_computed; + _system.eulerian_residual(jacobian_computed, context) && jacobian_computed; // Scale the time-dependent residual and jacobian correctly - _system.elem_residual *= _system.deltat; + context.elem_residual *= _system.deltat; if (jacobian_computed) - _system.elem_jacobian *= (theta * _system.deltat); + context.elem_jacobian *= (theta * _system.deltat); // The fixed_solution_derivative is always theta, // and now we're done scaling jacobians - _system.fixed_solution_derivative = theta; + context.fixed_solution_derivative = theta; // We evaluate mass_residual with the change in solution // to get the mass matrix, reusing old_elem_solution to hold that @@ -118,33 +119,33 @@ old_elem_solution -= theta_solution; // Move old_->elem_, theta_->old_ - _system.elem_solution.swap(old_elem_solution); + context.elem_solution.swap(old_elem_solution); // We do a trick here to avoid using a non-1 // elem_solution_derivative: - _system.elem_jacobian *= -1.0; - jacobian_computed = _system.mass_residual(jacobian_computed) && + context.elem_jacobian *= -1.0; + jacobian_computed = _system.mass_residual(jacobian_computed, context) && jacobian_computed; - _system.elem_jacobian *= -1.0; + context.elem_jacobian *= -1.0; // Move elem_->elem_, old_->theta_ - _system.elem_solution.swap(theta_solution); + context.elem_solution.swap(theta_solution); // Restore the elem position if necessary - _system.elem_reinit(1.); + context.elem_reinit(1.); // Add the constraint term - jacobian_computed = _system.element_constraint(jacobian_computed) && + jacobian_computed = _system.element_constraint(jacobian_computed, context) && jacobian_computed; // Add back the old residual and jacobian - _system.elem_residual += old_elem_residual; + context.elem_residual += old_elem_residual; if (request_jacobian) { if (jacobian_computed) - _system.elem_jacobian += old_elem_jacobian; + context.elem_jacobian += old_elem_jacobian; else - _system.elem_jacobian.swap(old_elem_jacobian); + context.elem_jacobian.swap(old_elem_jacobian); } return jacobian_computed; @@ -152,64 +153,65 @@ -bool EulerSolver::side_residual (bool request_jacobian) +bool EulerSolver::side_residual (bool request_jacobian, + DiffContext &context) { - unsigned int n_dofs = _system.elem_solution.size(); + unsigned int n_dofs = context.elem_solution.size(); // Local nonlinear solution at old timestep DenseVector old_elem_solution(n_dofs); for (unsigned int i=0; i != n_dofs; ++i) old_elem_solution(i) = - old_nonlinear_solution(_system.dof_indices[i]); + old_nonlinear_solution(context.dof_indices[i]); // Local nonlinear solution at time t_theta - DenseVector theta_solution(_system.elem_solution); + DenseVector theta_solution(context.elem_solution); theta_solution *= theta; theta_solution.add(1. - theta, old_elem_solution); // Technically the elem_solution_derivative is either theta // or 1.0 in this implementation, but we scale the former part // ourselves - _system.elem_solution_derivative = 1.0; + context.elem_solution_derivative = 1.0; // Technically the fixed_solution_derivative is always theta, // but we're scaling a whole jacobian by theta after these first // evaluations - _system.fixed_solution_derivative = 1.0; + context.fixed_solution_derivative = 1.0; // If a fixed solution is requested, we'll use theta_solution if (_system.use_fixed_solution) - _system.elem_fixed_solution = theta_solution; + context.elem_fixed_solution = theta_solution; // Move theta_->elem_, elem_->theta_ - _system.elem_solution.swap(theta_solution); + context.elem_solution.swap(theta_solution); // Move the mesh into place first if necessary - _system.elem_side_reinit(theta); + context.elem_side_reinit(theta); // We're going to compute just the change in elem_residual // (and possibly elem_jacobian), then add back the old values - DenseVector old_elem_residual(_system.elem_residual); + DenseVector old_elem_residual(context.elem_residual); DenseMatrix old_elem_jacobian; if (request_jacobian) { - old_elem_jacobian = _system.elem_jacobian; - _system.elem_jacobian.zero(); + old_elem_jacobian = context.elem_jacobian; + context.elem_jacobian.zero(); } - _system.elem_residual.zero(); + context.elem_residual.zero(); // Get the time derivative at t_theta bool jacobian_computed = - _system.side_time_derivative(request_jacobian); + _system.side_time_derivative(request_jacobian, context); // Scale the time-dependent residual and jacobian correctly - _system.elem_residual *= _system.deltat; + context.elem_residual *= _system.deltat; if (jacobian_computed) - _system.elem_jacobian *= (theta * _system.deltat); + context.elem_jacobian *= (theta * _system.deltat); // The fixed_solution_derivative is always theta, // and now we're done scaling jacobians - _system.fixed_solution_derivative = theta; + context.fixed_solution_derivative = theta; // We evaluate side_mass_residual with the change in solution // to get the mass matrix, reusing old_elem_solution to hold that @@ -219,33 +221,33 @@ old_elem_solution -= theta_solution; // Move old_->elem_, theta_->old_ - _system.elem_solution.swap(old_elem_solution); + context.elem_solution.swap(old_elem_solution); // We do a trick here to avoid using a non-1 // elem_solution_derivative: - _system.elem_jacobian *= -1.0; - jacobian_computed = _system.side_mass_residual(jacobian_computed) && + context.elem_jacobian *= -1.0; + jacobian_computed = _system.side_mass_residual(jacobian_computed, context) && jacobian_computed; - _system.elem_jacobian *= -1.0; + context.elem_jacobian *= -1.0; // Move elem_->elem_, old_->theta_ - _system.elem_solution.swap(theta_solution); + context.elem_solution.swap(theta_solution); // Restore the elem position if necessary - _system.elem_side_reinit(1.); + context.elem_side_reinit(1.); // Add the constraint term - jacobian_computed = _system.side_constraint(jacobian_computed) && + jacobian_computed = _system.side_constraint(jacobian_computed, context) && jacobian_computed; // Add back the old residual and jacobian - _system.elem_residual += old_elem_residual; + context.elem_residual += old_elem_residual; if (request_jacobian) { if (jacobian_computed) - _system.elem_jacobian += old_elem_jacobian; + context.elem_jacobian += old_elem_jacobian; else - _system.elem_jacobian.swap(old_elem_jacobian); + context.elem_jacobian.swap(old_elem_jacobian); } return jacobian_computed; Index: src/solvers/fem_context.C =================================================================== --- src/solvers/fem_context.C (revision 0) +++ src/solvers/fem_context.C (revision 0) @@ -0,0 +1,711 @@ +// $Id: fem_system.C 3118 2008-10-27 15:41:13Z jwpeterson $ + +// The libMesh Finite Element Library. +// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// 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 02111-1307 USA + + + +#include "dof_map.h" +#include "elem.h" +#include "fe_base.h" +#include "fe_interface.h" +#include "equation_systems.h" +#include "fem_context.h" +#include "fem_system.h" +#include "libmesh_logging.h" +#include "mesh.h" +#include "numeric_vector.h" +#include "quadrature.h" +#include "sparse_matrix.h" +#include "time_solver.h" +#include "unsteady_solver.h" // For euler_residual + + + + + +FEMContext::FEMContext (const FEMSystem &sys) + : DiffContext(sys), + element_qrule(NULL), side_qrule(NULL), + _mesh_sys(libMesh::invalid_uint), + _mesh_x_var(libMesh::invalid_uint), + _mesh_y_var(libMesh::invalid_uint), + _mesh_z_var(libMesh::invalid_uint), + elem(NULL), side(0), dim(sys.get_mesh().mesh_dimension()) +{ + // We need to know which of our variables has the hardest + // shape functions to numerically integrate. + + unsigned int n_vars = sys.n_vars(); + + libmesh_assert (n_vars); + FEType hardest_fe_type = sys.variable_type(0); + + for (unsigned int i=0; i != n_vars; ++i) + { + FEType fe_type = sys.variable_type(i); + + // FIXME - we don't yet handle mixed finite elements from + // different families which require different quadrature rules + // libmesh_assert (fe_type.family == hardest_fe_type.family); + + if (fe_type.order > hardest_fe_type.order) + hardest_fe_type = fe_type; + } + + // Create an adequate quadrature rule + element_qrule = hardest_fe_type.default_quadrature_rule + (dim, sys.extra_quadrature_order).release(); + side_qrule = hardest_fe_type.default_quadrature_rule + (dim-1, sys.extra_quadrature_order).release(); + + // Next, create finite element objects + element_fe_var.resize(n_vars); + side_fe_var.resize(n_vars); + for (unsigned int i=0; i != n_vars; ++i) + { + FEType fe_type = sys.variable_type(i); + if (element_fe[fe_type] == NULL) + { + element_fe[fe_type] = FEBase::build(dim, fe_type).release(); + element_fe[fe_type]->attach_quadrature_rule(element_qrule); + side_fe[fe_type] = FEBase::build(dim, fe_type).release(); + side_fe[fe_type]->attach_quadrature_rule(side_qrule); + } + element_fe_var[i] = element_fe[fe_type]; + side_fe_var[i] = side_fe[fe_type]; + } +} + + +FEMContext::~FEMContext() +{ + // We don't want to store AutoPtrs in STL containers, but we don't + // want to leak memory either + for (std::map::iterator i = element_fe.begin(); + i != element_fe.end(); ++i) + delete i->second; + element_fe.clear(); + + for (std::map::iterator i = side_fe.begin(); + i != side_fe.end(); ++i) + delete i->second; + side_fe.clear(); + + delete element_qrule; + element_qrule = NULL; + + delete side_qrule; + side_qrule = NULL; +} + + + +Number FEMContext::interior_value(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &phi = + element_fe_var[var]->get_phi(); + + // Accumulate solution value + Number u = 0.; + + for (unsigned int l=0; l != n_dofs; l++) + u += phi[l][qp] * coef(l); + + return u; +} + + + +Gradient FEMContext::interior_gradient(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &dphi = + element_fe_var[var]->get_dphi(); + + // Accumulate solution derivatives + Gradient du; + + for (unsigned int l=0; l != n_dofs; l++) + du.add_scaled(dphi[l][qp], coef(l)); + + return du; +} + + + +#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES +Tensor FEMContext::interior_hessian(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &d2phi = + element_fe_var[var]->get_d2phi(); + + // Accumulate solution second derivatives + Tensor d2u; + + for (unsigned int l=0; l != n_dofs; l++) + d2u.add_scaled(d2phi[l][qp], coef(l)); + + return d2u; +} +#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES + + + +Number FEMContext::side_value(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &phi = + side_fe_var[var]->get_phi(); + + // Accumulate solution value + Number u = 0.; + + for (unsigned int l=0; l != n_dofs; l++) + u += phi[l][qp] * coef(l); + + return u; +} + + + +Gradient FEMContext::side_gradient(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &dphi = + side_fe_var[var]->get_dphi(); + + // Accumulate solution derivatives + Gradient du; + + for (unsigned int l=0; l != n_dofs; l++) + du.add_scaled(dphi[l][qp], coef(l)); + + return du; +} + + + +#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES +Tensor FEMContext::side_hessian(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &d2phi = + side_fe_var[var]->get_d2phi(); + + // Accumulate solution second derivatives + Tensor d2u; + + for (unsigned int l=0; l != n_dofs; l++) + d2u.add_scaled(d2phi[l][qp], coef(l)); + + return d2u; +} +#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES + + + +Number FEMContext::point_value(unsigned int var, Point &p) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + + Number u = 0.; + + FEType fe_type = element_fe_var[var]->get_fe_type(); + Point p_master = FEInterface::inverse_map(dim, fe_type, elem, p); + + for (unsigned int l=0; l != n_dofs; l++) + u += FEInterface::shape(dim, fe_type, elem, l, p_master) + * coef(l); + + return u; +} + + + +Number FEMContext::fixed_interior_value(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_fixed_subsolutions.size() > var); + libmesh_assert (elem_fixed_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_fixed_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &phi = + element_fe_var[var]->get_phi(); + + // Accumulate solution value + Number u = 0.; + + for (unsigned int l=0; l != n_dofs; l++) + u += phi[l][qp] * coef(l); + + return u; +} + + + +Gradient FEMContext::fixed_interior_gradient(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_fixed_subsolutions.size() > var); + libmesh_assert (elem_fixed_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_fixed_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &dphi = + element_fe_var[var]->get_dphi(); + + // Accumulate solution derivatives + Gradient du; + + for (unsigned int l=0; l != n_dofs; l++) + du.add_scaled(dphi[l][qp], coef(l)); + + return du; +} + + + +#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES +Tensor FEMContext::fixed_interior_hessian(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_fixed_subsolutions.size() > var); + libmesh_assert (elem_fixed_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_fixed_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &d2phi = + element_fe_var[var]->get_d2phi(); + + // Accumulate solution second derivatives + Tensor d2u; + + for (unsigned int l=0; l != n_dofs; l++) + d2u.add_scaled(d2phi[l][qp], coef(l)); + + return d2u; +} +#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES + + + +Number FEMContext::fixed_side_value(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_fixed_subsolutions.size() > var); + libmesh_assert (elem_fixed_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_fixed_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &phi = + side_fe_var[var]->get_phi(); + + // Accumulate solution value + Number u = 0.; + + for (unsigned int l=0; l != n_dofs; l++) + u += phi[l][qp] * coef(l); + + return u; +} + + + +Gradient FEMContext::fixed_side_gradient(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_fixed_subsolutions.size() > var); + libmesh_assert (elem_fixed_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_fixed_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &dphi = + side_fe_var[var]->get_dphi(); + + // Accumulate solution derivatives + Gradient du; + + for (unsigned int l=0; l != n_dofs; l++) + du.add_scaled(dphi[l][qp], coef(l)); + + return du; +} + + + +#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES +Tensor FEMContext::fixed_side_hessian(unsigned int var, unsigned int qp) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + // Get current local coefficients + libmesh_assert (elem_fixed_subsolutions.size() > var); + libmesh_assert (elem_fixed_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_fixed_subsolutions[var]; + + // Get shape function values at quadrature point + const std::vector > &d2phi = + side_fe_var[var]->get_d2phi(); + + // Accumulate solution second derivatives + Tensor d2u; + + for (unsigned int l=0; l != n_dofs; l++) + d2u.add_scaled(d2phi[l][qp], coef(l)); + + return d2u; +} +#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES + + + +Number FEMContext::fixed_point_value(unsigned int var, Point &p) +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + + + // Get current local coefficients + libmesh_assert (elem_fixed_subsolutions.size() > var); + libmesh_assert (elem_fixed_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_fixed_subsolutions[var]; + + Number u = 0.; + + FEType fe_type = element_fe_var[var]->get_fe_type(); + Point p_master = FEInterface::inverse_map(dim, fe_type, elem, p); + + for (unsigned int l=0; l != n_dofs; l++) + u += FEInterface::shape(dim, fe_type, elem, l, p_master) + * coef(l); + + return u; +} + + + +void FEMContext::elem_reinit(Real theta) +{ + // Handle a moving element if necessary + if (_mesh_sys != libMesh::invalid_uint) + { + // FIXME - not threadsafe yet! + elem_position_set(theta); + elem_fe_reinit(); + } +} + + +void FEMContext::elem_side_reinit(Real theta) +{ + // Handle a moving element if necessary + if (_mesh_sys != libMesh::invalid_uint) + { + // FIXME - not threadsafe yet! + elem_position_set(theta); + elem_side_fe_reinit(); + } +} + + +void FEMContext::elem_fe_reinit () +{ + // Initialize all the interior FE objects on elem. + // Logging of FE::reinit is done in the FE functions + std::map::iterator fe_end = element_fe.end(); + for (std::map::iterator i = element_fe.begin(); + i != fe_end; ++i) + { + i->second->reinit(elem); + } +} + + +void FEMContext::elem_side_fe_reinit () +{ + // Initialize all the interior FE objects on elem/side. + // Logging of FE::reinit is done in the FE functions + std::map::iterator fe_end = side_fe.end(); + for (std::map::iterator i = side_fe.begin(); + i != fe_end; ++i) + { + i->second->reinit(elem, side); + } +} + + + +void FEMContext::elem_position_get() +{ + // This is too expensive to call unless we've been asked to move the mesh + libmesh_assert (_mesh_sys != libMesh::invalid_uint); + + // This will probably break with threading when two contexts are + // operating on elements which share a node + untested(); + + // If the coordinate data is in our own system, it's already + // been set up for us +// if (_mesh_sys == this->number()) +// { + unsigned int n_nodes = elem->n_nodes(); + // For simplicity we demand that mesh coordinates be stored + // in a format that allows a direct copy + libmesh_assert(_mesh_x_var == libMesh::invalid_uint || + (element_fe_var[_mesh_x_var]->get_fe_type().family + == LAGRANGE && + element_fe_var[_mesh_x_var]->get_fe_type().order + == elem->default_order())); + libmesh_assert(_mesh_y_var == libMesh::invalid_uint || + (element_fe_var[_mesh_y_var]->get_fe_type().family + == LAGRANGE && + element_fe_var[_mesh_y_var]->get_fe_type().order + == elem->default_order())); + libmesh_assert(_mesh_z_var == libMesh::invalid_uint || + (element_fe_var[_mesh_z_var]->get_fe_type().family + == LAGRANGE && + element_fe_var[_mesh_z_var]->get_fe_type().order + == elem->default_order())); + + // Get degree of freedom coefficients from point coordinates + if (_mesh_x_var != libMesh::invalid_uint) + for (unsigned int i=0; i != n_nodes; ++i) + (*elem_subsolutions[_mesh_x_var])(i) = elem->point(i)(0); + + if (_mesh_y_var != libMesh::invalid_uint) + for (unsigned int i=0; i != n_nodes; ++i) + (*elem_subsolutions[_mesh_y_var])(i) = elem->point(i)(1); + + if (_mesh_z_var != libMesh::invalid_uint) + for (unsigned int i=0; i != n_nodes; ++i) + (*elem_subsolutions[_mesh_z_var])(i) = elem->point(i)(2); +// } + // FIXME - If the coordinate data is not in our own system, someone + // had better get around to implementing that... - RHS +// else +// { +// libmesh_not_implemented(); +// } +} + + + +void FEMContext::elem_position_set(Real) +{ + // This is too expensive to call unless we've been asked to move the mesh + libmesh_assert (_mesh_sys != libMesh::invalid_uint); + + // This will probably break with threading when two contexts are + // operating on elements which share a node + untested(); + + // If the coordinate data is in our own system, it's already + // been set up for us, and we can ignore our input parameter theta +// if (_mesh_sys == this->number()) +// { + unsigned int n_nodes = elem->n_nodes(); + // For simplicity we demand that mesh coordinates be stored + // in a format that allows a direct copy + libmesh_assert(_mesh_x_var == libMesh::invalid_uint || + (element_fe_var[_mesh_x_var]->get_fe_type().family + == LAGRANGE && + elem_subsolutions[_mesh_x_var]->size() == n_nodes)); + libmesh_assert(_mesh_y_var == libMesh::invalid_uint || + (element_fe_var[_mesh_y_var]->get_fe_type().family + == LAGRANGE && + elem_subsolutions[_mesh_y_var]->size() == n_nodes)); + libmesh_assert(_mesh_z_var == libMesh::invalid_uint || + (element_fe_var[_mesh_z_var]->get_fe_type().family + == LAGRANGE && + elem_subsolutions[_mesh_z_var]->size() == n_nodes)); + + // Set the new point coordinates + if (_mesh_x_var != libMesh::invalid_uint) + for (unsigned int i=0; i != n_nodes; ++i) + elem->point(i)(0) = + libmesh_real((*elem_subsolutions[_mesh_x_var])(i)); + + if (_mesh_y_var != libMesh::invalid_uint) + for (unsigned int i=0; i != n_nodes; ++i) + elem->point(i)(1) = + libmesh_real((*elem_subsolutions[_mesh_y_var])(i)); + + if (_mesh_z_var != libMesh::invalid_uint) + for (unsigned int i=0; i != n_nodes; ++i) + elem->point(i)(2) = + libmesh_real((*elem_subsolutions[_mesh_z_var])(i)); +// } + // FIXME - If the coordinate data is not in our own system, someone + // had better get around to implementing that... - RHS +// else +// { +// libmesh_not_implemented(); +// } +} + + + + + +void FEMContext::reinit(const FEMSystem &sys, Elem *e) +{ + elem = e; + // Initialize the per-element data for elem. + sys.get_dof_map().dof_indices (elem, dof_indices); + unsigned int n_dofs = dof_indices.size(); + + elem_solution.resize(n_dofs); + if (sys.use_fixed_solution) + elem_fixed_solution.resize(n_dofs); + + for (unsigned int i=0; i != n_dofs; ++i) + elem_solution(i) = sys.current_solution(dof_indices[i]); + + // These resize calls also zero out the residual and jacobian + elem_residual.resize(n_dofs); + elem_jacobian.resize(n_dofs, n_dofs); + + // Initialize the per-variable data for elem. + unsigned int sub_dofs = 0; + for (unsigned int i=0; i != sys.n_vars(); ++i) + { + sys.get_dof_map().dof_indices (elem, dof_indices_var[i], i); + + elem_subsolutions[i]->reposition + (sub_dofs, dof_indices_var[i].size()); + + if (sys.use_fixed_solution) + elem_fixed_subsolutions[i]->reposition + (sub_dofs, dof_indices_var[i].size()); + + elem_subresiduals[i]->reposition + (sub_dofs, dof_indices_var[i].size()); + + for (unsigned int j=0; j != i; ++j) + { + elem_subjacobians[i][j]->reposition + (sub_dofs, elem_subresiduals[j]->i_off(), + dof_indices_var[i].size(), + dof_indices_var[j].size()); + elem_subjacobians[j][i]->reposition + (elem_subresiduals[j]->i_off(), sub_dofs, + dof_indices_var[j].size(), + dof_indices_var[i].size()); + } + elem_subjacobians[i][i]->reposition + (sub_dofs, sub_dofs, + dof_indices_var[i].size(), + dof_indices_var[i].size()); + sub_dofs += dof_indices_var[i].size(); + } + libmesh_assert(sub_dofs == n_dofs); + + // Moving the mesh may be necessary + if (_mesh_sys != libMesh::invalid_uint) + elem_position_set(1.); + // Reinitializing the FE objects is definitely necessary + this->elem_fe_reinit(); +} Index: src/solvers/fem_system.C =================================================================== --- src/solvers/fem_system.C (revision 3205) +++ src/solvers/fem_system.C (working copy) @@ -45,7 +45,6 @@ extra_quadrature_order(0), numerical_jacobian_h(1.e-6), verify_analytic_jacobians(0.0), - element_qrule(NULL), side_qrule(NULL), elem(NULL), _mesh_sys(libMesh::invalid_uint), _mesh_x_var(libMesh::invalid_uint), _mesh_y_var(libMesh::invalid_uint), @@ -61,517 +60,20 @@ -Number FEMSystem::interior_value(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_subsolutions.size() > var); - libmesh_assert (elem_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &phi = - element_fe_var[var]->get_phi(); - - // Accumulate solution value - Number u = 0.; - - for (unsigned int l=0; l != n_dofs; l++) - u += phi[l][qp] * coef(l); - - return u; -} - - - -Gradient FEMSystem::interior_gradient(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_subsolutions.size() > var); - libmesh_assert (elem_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &dphi = - element_fe_var[var]->get_dphi(); - - // Accumulate solution derivatives - Gradient du; - - for (unsigned int l=0; l != n_dofs; l++) - du.add_scaled(dphi[l][qp], coef(l)); - - return du; -} - - - -#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES -Tensor FEMSystem::interior_hessian(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_subsolutions.size() > var); - libmesh_assert (elem_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &d2phi = - element_fe_var[var]->get_d2phi(); - - // Accumulate solution second derivatives - Tensor d2u; - - for (unsigned int l=0; l != n_dofs; l++) - d2u.add_scaled(d2phi[l][qp], coef(l)); - - return d2u; -} -#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES - - - -Number FEMSystem::side_value(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_subsolutions.size() > var); - libmesh_assert (elem_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &phi = - side_fe_var[var]->get_phi(); - - // Accumulate solution value - Number u = 0.; - - for (unsigned int l=0; l != n_dofs; l++) - u += phi[l][qp] * coef(l); - - return u; -} - - - -Gradient FEMSystem::side_gradient(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_subsolutions.size() > var); - libmesh_assert (elem_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &dphi = - side_fe_var[var]->get_dphi(); - - // Accumulate solution derivatives - Gradient du; - - for (unsigned int l=0; l != n_dofs; l++) - du.add_scaled(dphi[l][qp], coef(l)); - - return du; -} - - - -#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES -Tensor FEMSystem::side_hessian(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_subsolutions.size() > var); - libmesh_assert (elem_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &d2phi = - side_fe_var[var]->get_d2phi(); - - // Accumulate solution second derivatives - Tensor d2u; - - for (unsigned int l=0; l != n_dofs; l++) - d2u.add_scaled(d2phi[l][qp], coef(l)); - - return d2u; -} -#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES - - - -Number FEMSystem::point_value(unsigned int var, Point &p) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - - // Get current local coefficients - libmesh_assert (elem_subsolutions.size() > var); - libmesh_assert (elem_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_subsolutions[var]; - - Number u = 0.; - - unsigned int dim = get_mesh().mesh_dimension(); - FEType fe_type = element_fe_var[var]->get_fe_type(); - Point p_master = FEInterface::inverse_map(dim, fe_type, elem, p); - - for (unsigned int l=0; l != n_dofs; l++) - u += FEInterface::shape(dim, fe_type, elem, l, p_master) - * coef(l); - - return u; -} - - - -Number FEMSystem::fixed_interior_value(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_fixed_subsolutions.size() > var); - libmesh_assert (elem_fixed_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_fixed_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &phi = - element_fe_var[var]->get_phi(); - - // Accumulate solution value - Number u = 0.; - - for (unsigned int l=0; l != n_dofs; l++) - u += phi[l][qp] * coef(l); - - return u; -} - - - -Gradient FEMSystem::fixed_interior_gradient(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_fixed_subsolutions.size() > var); - libmesh_assert (elem_fixed_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_fixed_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &dphi = - element_fe_var[var]->get_dphi(); - - // Accumulate solution derivatives - Gradient du; - - for (unsigned int l=0; l != n_dofs; l++) - du.add_scaled(dphi[l][qp], coef(l)); - - return du; -} - - - -#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES -Tensor FEMSystem::fixed_interior_hessian(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_fixed_subsolutions.size() > var); - libmesh_assert (elem_fixed_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_fixed_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &d2phi = - element_fe_var[var]->get_d2phi(); - - // Accumulate solution second derivatives - Tensor d2u; - - for (unsigned int l=0; l != n_dofs; l++) - d2u.add_scaled(d2phi[l][qp], coef(l)); - - return d2u; -} -#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES - - - -Number FEMSystem::fixed_side_value(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_fixed_subsolutions.size() > var); - libmesh_assert (elem_fixed_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_fixed_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &phi = - side_fe_var[var]->get_phi(); - - // Accumulate solution value - Number u = 0.; - - for (unsigned int l=0; l != n_dofs; l++) - u += phi[l][qp] * coef(l); - - return u; -} - - - -Gradient FEMSystem::fixed_side_gradient(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_fixed_subsolutions.size() > var); - libmesh_assert (elem_fixed_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_fixed_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &dphi = - side_fe_var[var]->get_dphi(); - - // Accumulate solution derivatives - Gradient du; - - for (unsigned int l=0; l != n_dofs; l++) - du.add_scaled(dphi[l][qp], coef(l)); - - return du; -} - - - -#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES -Tensor FEMSystem::fixed_side_hessian(unsigned int var, unsigned int qp) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - // Get current local coefficients - libmesh_assert (elem_fixed_subsolutions.size() > var); - libmesh_assert (elem_fixed_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_fixed_subsolutions[var]; - - // Get shape function values at quadrature point - const std::vector > &d2phi = - side_fe_var[var]->get_d2phi(); - - // Accumulate solution second derivatives - Tensor d2u; - - for (unsigned int l=0; l != n_dofs; l++) - d2u.add_scaled(d2phi[l][qp], coef(l)); - - return d2u; -} -#endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES - - - -Number FEMSystem::fixed_point_value(unsigned int var, Point &p) -{ - // Get local-to-global dof index lookup - libmesh_assert (dof_indices.size() > var); - const unsigned int n_dofs = dof_indices_var[var].size(); - - - // Get current local coefficients - libmesh_assert (elem_fixed_subsolutions.size() > var); - libmesh_assert (elem_fixed_subsolutions[var] != NULL); - DenseSubVector &coef = *elem_fixed_subsolutions[var]; - - Number u = 0.; - - unsigned int dim = get_mesh().mesh_dimension(); - FEType fe_type = element_fe_var[var]->get_fe_type(); - Point p_master = FEInterface::inverse_map(dim, fe_type, elem, p); - - for (unsigned int l=0; l != n_dofs; l++) - u += FEInterface::shape(dim, fe_type, elem, l, p_master) - * coef(l); - - return u; -} - - - void FEMSystem::clear() { - this->clear_fem_ptrs(); - Parent::clear(); } -void FEMSystem::clear_fem_ptrs() -{ - // We don't want to store AutoPtrs in STL containers, but we don't - // want to leak memory either - for (std::map::iterator i = element_fe.begin(); - i != element_fe.end(); ++i) - delete i->second; - element_fe.clear(); - - for (std::map::iterator i = side_fe.begin(); - i != side_fe.end(); ++i) - delete i->second; - side_fe.clear(); - - delete element_qrule; - element_qrule = NULL; - - delete side_qrule; - side_qrule = NULL; -} - - - void FEMSystem::init_data () { - // We may have already been initialized once - this->clear_fem_ptrs(); - // First initialize LinearImplicitSystem data Parent::init_data(); - - const MeshBase &mesh = this->get_mesh(); - - unsigned int dim = mesh.mesh_dimension(); - - // We need to know which of our variables has the hardest - // shape functions to numerically integrate. - - unsigned int n_vars = this->n_vars(); - - libmesh_assert (n_vars); - FEType hardest_fe_type = this->variable_type(0); - - for (unsigned int i=0; i != n_vars; ++i) - { - FEType fe_type = this->variable_type(i); - - // FIXME - we don't yet handle mixed finite elements from - // different families which require different quadrature rules - // libmesh_assert (fe_type.family == hardest_fe_type.family); - - if (fe_type.order > hardest_fe_type.order) - hardest_fe_type = fe_type; - } - - // Create an adequate quadrature rule - element_qrule = hardest_fe_type.default_quadrature_rule - (dim, extra_quadrature_order).release(); - side_qrule = hardest_fe_type.default_quadrature_rule - (dim-1, extra_quadrature_order).release(); - - // Next, create finite element objects - element_fe_var.resize(n_vars); - side_fe_var.resize(n_vars); - for (unsigned int i=0; i != n_vars; ++i) - { - FEType fe_type = this->variable_type(i); - if (element_fe[fe_type] == NULL) - { - element_fe[fe_type] = FEBase::build(dim, fe_type).release(); - element_fe[fe_type]->attach_quadrature_rule(element_qrule); - side_fe[fe_type] = FEBase::build(dim, fe_type).release(); - side_fe[fe_type]->attach_quadrature_rule(side_qrule); - } - element_fe_var[i] = element_fe[fe_type]; - side_fe_var[i] = side_fe[fe_type]; - } } -void FEMSystem::elem_reinit(Real theta) -{ - // Handle a moving element if necessary - if (_mesh_sys != libMesh::invalid_uint) - { - elem_position_set(theta); - elem_fe_reinit(); - } -} - - -void FEMSystem::elem_side_reinit(Real theta) -{ - // Handle a moving element if necessary - if (_mesh_sys != libMesh::invalid_uint) - { - elem_position_set(theta); - elem_side_fe_reinit(); - } -} - - -void FEMSystem::elem_fe_reinit () -{ - // Initialize all the interior FE objects on elem. - // Logging of FE::reinit is done in the FE functions - std::map::iterator fe_end = element_fe.end(); - for (std::map::iterator i = element_fe.begin(); - i != fe_end; ++i) - { - i->second->reinit(elem); - } -} - - -void FEMSystem::elem_side_fe_reinit () -{ - // Initialize all the interior FE objects on elem/side. - // Logging of FE::reinit is done in the FE functions - std::map::iterator fe_end = side_fe.end(); - for (std::map::iterator i = side_fe.begin(); - i != fe_end; ++i) - { - i->second->reinit(elem, side); - } -} - - void FEMSystem::assembly (bool get_residual, bool get_jacobian) { libmesh_assert(get_residual || get_jacobian); @@ -587,8 +89,6 @@ const MeshBase& mesh = this->get_mesh(); - const DofMap& dof_map = this->get_dof_map(); - // this->get_vector("_nonlinear_solution").localize // (*current_local_nonlinear_solution, // dof_map.get_send_list()); @@ -634,6 +134,9 @@ // we're using libmesh_assert (time_solver.get() != NULL); + AutoPtr con = this->build_context(); + FEMContext &_femcontext = libmesh_cast_ref(*con); + // Build the residual and jacobian contributions on every active // mesh element on this processor MeshBase::const_element_iterator el = @@ -643,74 +146,18 @@ for ( ; el != end_el; ++el) { - elem = *el; + _femcontext.reinit(*this, *el); - // Initialize the per-element data for elem. - dof_map.dof_indices (elem, dof_indices); - unsigned int n_dofs = dof_indices.size(); + bool jacobian_computed = + time_solver->element_residual(get_jacobian, _femcontext); - elem_solution.resize(n_dofs); - if (use_fixed_solution) - elem_fixed_solution.resize(n_dofs); - - for (unsigned int i=0; i != n_dofs; ++i) -// elem_solution(i) = current_nonlinear_solution(dof_indices[i]); - elem_solution(i) = current_solution(dof_indices[i]); - - // These resize calls also zero out the residual and jacobian - elem_residual.resize(n_dofs); - elem_jacobian.resize(n_dofs, n_dofs); - - // Initialize the per-variable data for elem. - unsigned int sub_dofs = 0; - for (unsigned int i=0; i != this->n_vars(); ++i) - { - dof_map.dof_indices (elem, dof_indices_var[i], i); - - elem_subsolutions[i]->reposition - (sub_dofs, dof_indices_var[i].size()); - - if (use_fixed_solution) - elem_fixed_subsolutions[i]->reposition - (sub_dofs, dof_indices_var[i].size()); - - elem_subresiduals[i]->reposition - (sub_dofs, dof_indices_var[i].size()); - - for (unsigned int j=0; j != i; ++j) - { - elem_subjacobians[i][j]->reposition - (sub_dofs, elem_subresiduals[j]->i_off(), - dof_indices_var[i].size(), - dof_indices_var[j].size()); - elem_subjacobians[j][i]->reposition - (elem_subresiduals[j]->i_off(), sub_dofs, - dof_indices_var[j].size(), - dof_indices_var[i].size()); - } - elem_subjacobians[i][i]->reposition - (sub_dofs, sub_dofs, - dof_indices_var[i].size(), - dof_indices_var[i].size()); - sub_dofs += dof_indices_var[i].size(); - } - libmesh_assert(sub_dofs == n_dofs); - - // Moving the mesh may be necessary - if (_mesh_sys != libMesh::invalid_uint) - elem_position_set(1.); - // Reinitializing the FE objects is definitely necessary - this->elem_fe_reinit(); - - bool jacobian_computed = time_solver->element_residual(get_jacobian); - // Compute a numeric jacobian if we have to if (get_jacobian && !jacobian_computed) { // Make sure we didn't compute a jacobian and lie about it - libmesh_assert(elem_jacobian.l1_norm() == 0.0); + libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0); // Logging of numerical jacobians is done separately - this->numerical_elem_jacobian(); + this->numerical_elem_jacobian(_femcontext); } // Compute a numeric jacobian if we're asked to verify the @@ -718,20 +165,20 @@ if (get_jacobian && jacobian_computed && verify_analytic_jacobians != 0.0) { - DenseMatrix analytic_jacobian(elem_jacobian); + DenseMatrix analytic_jacobian(_femcontext.elem_jacobian); - elem_jacobian.zero(); + _femcontext.elem_jacobian.zero(); // Logging of numerical jacobians is done separately - this->numerical_elem_jacobian(); + this->numerical_elem_jacobian(_femcontext); Real analytic_norm = analytic_jacobian.l1_norm(); - Real numerical_norm = elem_jacobian.l1_norm(); + Real numerical_norm = _femcontext.elem_jacobian.l1_norm(); // If we can continue, we'll probably prefer the analytic jacobian - analytic_jacobian.swap(elem_jacobian); + analytic_jacobian.swap(_femcontext.elem_jacobian); // The matrix "analytic_jacobian" will now hold the error matrix - analytic_jacobian.add(-1.0, elem_jacobian); + analytic_jacobian.add(-1.0, _femcontext.elem_jacobian); Real error_norm = analytic_jacobian.l1_norm(); Real relative_error = error_norm / @@ -741,14 +188,14 @@ { std::cerr << "Relative error " << relative_error << " detected in analytic jacobian on element " - << elem->id() << '!' << std::endl; + << _femcontext.elem->id() << '!' << std::endl; unsigned int old_precision = std::cout.precision(); std::cout.precision(16); - std::cout << "J_analytic " << elem->id() << " = " - << elem_jacobian << std::endl; - analytic_jacobian.add(1.0, elem_jacobian); - std::cout << "J_numeric " << elem->id() << " = " + std::cout << "J_analytic " << _femcontext.elem->id() << " = " + << _femcontext.elem_jacobian << std::endl; + analytic_jacobian.add(1.0, _femcontext.elem_jacobian); + std::cout << "J_numeric " << _femcontext.elem->id() << " = " << analytic_jacobian << std::endl; std::cout.precision(old_precision); @@ -757,16 +204,19 @@ } } - for (side = 0; side != elem->n_sides(); ++side) + for (_femcontext.side = 0; + _femcontext.side != _femcontext.elem->n_sides(); + ++_femcontext.side) { // Don't compute on non-boundary sides unless requested - if (!compute_internal_sides && elem->neighbor(side) != NULL) + if (!compute_internal_sides && + _femcontext.elem->neighbor(_femcontext.side) != NULL) continue; // Any mesh movement has already been done (and restored, // if the TimeSolver isn't broken), but // reinitializing the side FE objects is still necessary - this->elem_side_fe_reinit(); + _femcontext.elem_side_fe_reinit(); DenseMatrix old_jacobian; // If we're in DEBUG mode, we should always verify that the @@ -779,10 +229,11 @@ if (verify_analytic_jacobians != 0.0 && get_jacobian) #endif // ifndef DEBUG { - old_jacobian = elem_jacobian; - elem_jacobian.zero(); + old_jacobian = _femcontext.elem_jacobian; + _femcontext.elem_jacobian.zero(); } - jacobian_computed = time_solver->side_residual(get_jacobian); + jacobian_computed = + time_solver->side_residual(get_jacobian, _femcontext); // Compute a numeric jacobian if we have to if (get_jacobian && !jacobian_computed) @@ -791,10 +242,10 @@ // so we can make sure side_residual didn't compute a // jacobian and lie about it #ifdef DEBUG - libmesh_assert(elem_jacobian.l1_norm() == 0.0); + libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0); #endif // Logging of numerical jacobians is done separately - this->numerical_side_jacobian(); + this->numerical_side_jacobian(_femcontext); // If we're in DEBUG mode or if // verify_analytic_jacobians is on, we've moved @@ -803,7 +254,7 @@ #ifndef DEBUG if (verify_analytic_jacobians != 0.0) #endif // ifndef DEBUG - elem_jacobian += old_jacobian; + _femcontext.elem_jacobian += old_jacobian; } // Compute a numeric jacobian if we're asked to verify the @@ -811,20 +262,20 @@ if (get_jacobian && jacobian_computed && verify_analytic_jacobians != 0.0) { - DenseMatrix analytic_jacobian(elem_jacobian); + DenseMatrix analytic_jacobian(_femcontext.elem_jacobian); - elem_jacobian.zero(); + _femcontext.elem_jacobian.zero(); // Logging of numerical jacobians is done separately - this->numerical_side_jacobian(); + this->numerical_side_jacobian(_femcontext); Real analytic_norm = analytic_jacobian.l1_norm(); - Real numerical_norm = elem_jacobian.l1_norm(); + Real numerical_norm = _femcontext.elem_jacobian.l1_norm(); // If we can continue, we'll probably prefer the analytic jacobian - analytic_jacobian.swap(elem_jacobian); + analytic_jacobian.swap(_femcontext.elem_jacobian); // The matrix "analytic_jacobian" will now hold the error matrix - analytic_jacobian.add(-1.0, elem_jacobian); + analytic_jacobian.add(-1.0, _femcontext.elem_jacobian); Real error_norm = analytic_jacobian.l1_norm(); Real relative_error = error_norm / @@ -834,15 +285,16 @@ { std::cerr << "Relative error " << relative_error << " detected in analytic jacobian on element " - << elem->id() - << ", side " << side << '!' << std::endl; + << _femcontext.elem->id() + << ", side " + << _femcontext.side << '!' << std::endl; unsigned int old_precision = std::cout.precision(); std::cout.precision(16); - std::cout << "J_analytic " << elem->id() << " = " - << elem_jacobian << std::endl; - analytic_jacobian.add(1.0, elem_jacobian); - std::cout << "J_numeric " << elem->id() << " = " + std::cout << "J_analytic " << _femcontext.elem->id() << " = " + << _femcontext.elem_jacobian << std::endl; + analytic_jacobian.add(1.0, _femcontext.elem_jacobian); + std::cout << "J_numeric " << _femcontext.elem->id() << " = " << analytic_jacobian << std::endl; std::cout.precision(old_precision); @@ -850,7 +302,7 @@ } // Once we've verified a side, we'll want to add back the // rest of the accumulated jacobian - elem_jacobian += old_jacobian; + _femcontext.elem_jacobian += old_jacobian; } // In DEBUG mode, we've set elem_jacobian == 0, and we // may still need to add the old jacobian back @@ -858,7 +310,7 @@ if (get_jacobian && jacobian_computed && verify_analytic_jacobians == 0.0) { - elem_jacobian += old_jacobian; + _femcontext.elem_jacobian += old_jacobian; } #endif // ifdef DEBUG } @@ -868,28 +320,31 @@ // enforce_constraints_exactly() should be called in the solver if (get_residual && get_jacobian) this->get_dof_map().constrain_element_matrix_and_vector - (elem_jacobian, elem_residual, dof_indices, false); + (_femcontext.elem_jacobian, _femcontext.elem_residual, + _femcontext.dof_indices, false); else if (get_residual) this->get_dof_map().constrain_element_vector - (elem_residual, dof_indices, false); + (_femcontext.elem_residual, _femcontext.dof_indices, false); else if (get_jacobian) this->get_dof_map().constrain_element_matrix - (elem_jacobian, dof_indices, false); + (_femcontext.elem_jacobian, _femcontext.dof_indices, false); #endif // #ifdef LIBMESH_ENABLE_AMR if (get_jacobian && print_element_jacobians) { unsigned int old_precision = std::cout.precision(); std::cout.precision(16); - std::cout << "J_elem " << elem->id() << " = " - << elem_jacobian << std::endl; + std::cout << "J_elem " << _femcontext.elem->id() << " = " + << _femcontext.elem_jacobian << std::endl; std::cout.precision(old_precision); } if (get_jacobian) - this->matrix->add_matrix (elem_jacobian, dof_indices); + this->matrix->add_matrix (_femcontext.elem_jacobian, + _femcontext.dof_indices); if (get_residual) - this->rhs->add_vector (elem_residual, dof_indices); + this->rhs->add_vector (_femcontext.elem_residual, + _femcontext.dof_indices); } @@ -941,12 +396,13 @@ const MeshBase& mesh = this->get_mesh(); - const DofMap& dof_map = this->get_dof_map(); - // This code won't work on a parallelized mesh yet - // it won't get ancestor elements right. libmesh_assert(mesh.is_serial()); + AutoPtr con = this->build_context(); + FEMContext &_femcontext = libmesh_cast_ref(*con); + // Move every mesh element we can MeshBase::const_element_iterator el = mesh.active_elements_begin(); @@ -955,33 +411,12 @@ for ( ; el != end_el; ++el) { - elem = *el; + _femcontext.reinit(*this, *el); // This code won't handle moving subactive elements - libmesh_assert(!elem->has_children()); + libmesh_assert(!_femcontext.elem->has_children()); - // Initialize the per-element data for elem. - dof_map.dof_indices (elem, dof_indices); - unsigned int n_dofs = dof_indices.size(); - - elem_solution.resize(n_dofs); - - for (unsigned int i=0; i != n_dofs; ++i) - elem_solution(i) = current_solution(dof_indices[i]); - - // Initialize the per-variable data for elem. - unsigned int sub_dofs = 0; - for (unsigned int i=0; i != this->n_vars(); ++i) - { - dof_map.dof_indices (elem, dof_indices_var[i], i); - - elem_subsolutions[i]->reposition - (sub_dofs, dof_indices_var[i].size()); - - sub_dofs += dof_indices_var[i].size(); - } - - this->elem_position_set(0.); + _femcontext.elem_position_set(0.); } } @@ -993,13 +428,11 @@ const MeshBase& mesh = this->get_mesh(); - const DofMap& dof_map = this->get_dof_map(); - -// this->get_vector("_nonlinear_solution").localize -// (*current_local_nonlinear_solution, -// dof_map.get_send_list()); this->update(); + AutoPtr con = this->build_context(); + FEMContext &_femcontext = libmesh_cast_ref(*con); + // Loop over every active mesh element on this processor MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); @@ -1008,68 +441,40 @@ for ( ; el != end_el; ++el) { - elem = *el; + _femcontext.reinit(*this, *el); - // Initialize the per-element data for elem. - dof_map.dof_indices (elem, dof_indices); - unsigned int n_dofs = dof_indices.size(); - - elem_solution.resize(n_dofs); - // This resize call also zeros out the residual - elem_residual.resize(n_dofs); - - for (unsigned int i=0; i != n_dofs; ++i) -// elem_solution(i) = current_nonlinear_solution(dof_indices[i]); - elem_solution(i) = current_solution(dof_indices[i]); - - // Initialize the per-variable data for elem. - unsigned int sub_dofs = 0; - for (unsigned int i=0; i != this->n_vars(); ++i) - { - dof_map.dof_indices (elem, dof_indices_var[i], i); - - elem_subsolutions[i]->reposition - (sub_dofs, dof_indices_var[i].size()); - - if (use_fixed_solution) - elem_fixed_subsolutions[i]->reposition - (sub_dofs, dof_indices_var[i].size()); - - elem_subresiduals[i]->reposition - (sub_dofs, dof_indices_var[i].size()); - - sub_dofs += dof_indices_var[i].size(); - } - libmesh_assert(sub_dofs == n_dofs); - // Optionally initialize all the interior FE objects on elem. - if (fe_reinit_during_postprocess) - this->elem_fe_reinit(); + // if (fe_reinit_during_postprocess) + // _femcontext.elem_fe_reinit(); - this->element_postprocess(); + this->element_postprocess(_femcontext); - for (side = 0; side != elem->n_sides(); ++side) + for (_femcontext.side = 0; + _femcontext.side != _femcontext.elem->n_sides(); + ++_femcontext.side) { // Don't compute on non-boundary sides unless requested if (!postprocess_sides || (!compute_internal_sides && - elem->neighbor(side) != NULL)) + _femcontext.elem->neighbor(_femcontext.side) != NULL)) continue; // Optionally initialize all the interior FE objects on elem/side. // Logging of FE::reinit is done in the FE functions if (fe_reinit_during_postprocess) { - std::map::iterator fe_end = element_fe.end(); - fe_end = side_fe.end(); - for (std::map::iterator i = side_fe.begin(); + std::map::iterator fe_end = + _femcontext.element_fe.end(); + fe_end = _femcontext.side_fe.end(); + for (std::map::iterator i = + _femcontext.side_fe.begin(); i != fe_end; ++i) { - i->second->reinit(elem, side); + i->second->reinit(_femcontext.elem, _femcontext.side); } } - this->side_postprocess(); + this->side_postprocess(_femcontext); } } @@ -1078,27 +483,28 @@ -void FEMSystem::numerical_jacobian (TimeSolverResPtr res) +void FEMSystem::numerical_jacobian (TimeSolverResPtr res, + FEMContext &context) { // Logging is done by numerical_elem_jacobian // or numerical_side_jacobian - DenseVector original_residual(elem_residual); - DenseVector backwards_residual(elem_residual); - DenseMatrix numerical_jacobian(elem_jacobian); + DenseVector original_residual(context.elem_residual); + DenseVector backwards_residual(context.elem_residual); + DenseMatrix numerical_jacobian(context.elem_jacobian); #ifdef DEBUG - DenseMatrix old_jacobian(elem_jacobian); + DenseMatrix old_jacobian(context.elem_jacobian); #endif Real numerical_point_h = 0.; if (_mesh_sys == this->number()) - numerical_point_h = numerical_jacobian_h * elem->hmin(); + numerical_point_h = numerical_jacobian_h * context.elem->hmin(); - for (unsigned int j = 0; j != dof_indices.size(); ++j) + for (unsigned int j = 0; j != context.dof_indices.size(); ++j) { // Take the "minus" side of a central differenced first derivative - Number original_solution = elem_solution(j); - elem_solution(j) -= numerical_jacobian_h; + Number original_solution = context.elem_solution(j); + context.elem_solution(j) -= numerical_jacobian_h; // Make sure to catch any moving mesh terms // FIXME - this could be less ugly @@ -1107,102 +513,125 @@ { if (_mesh_x_var != libMesh::invalid_uint) for (unsigned int k = 0; - k != dof_indices_var[_mesh_x_var].size(); ++k) - if (dof_indices_var[_mesh_x_var][k] == dof_indices[j]) - coord = &(elem->point(k)(0)); + k != context.dof_indices_var[_mesh_x_var].size(); ++k) + if (context.dof_indices_var[_mesh_x_var][k] == + context.dof_indices[j]) + coord = &(context.elem->point(k)(0)); if (_mesh_y_var != libMesh::invalid_uint) for (unsigned int k = 0; - k != dof_indices_var[_mesh_y_var].size(); ++k) - if (dof_indices_var[_mesh_y_var][k] == dof_indices[j]) - coord = &(elem->point(k)(1)); + k != context.dof_indices_var[_mesh_y_var].size(); ++k) + if (context.dof_indices_var[_mesh_y_var][k] == + context.dof_indices[j]) + coord = &(context.elem->point(k)(1)); if (_mesh_z_var != libMesh::invalid_uint) for (unsigned int k = 0; - k != dof_indices_var[_mesh_z_var].size(); ++k) - if (dof_indices_var[_mesh_z_var][k] == dof_indices[j]) - coord = &(elem->point(k)(2)); + k != context.dof_indices_var[_mesh_z_var].size(); ++k) + if (context.dof_indices_var[_mesh_z_var][k] == + context.dof_indices[j]) + coord = &(context.elem->point(k)(2)); } if (coord) { // We have enough information to scale the perturbations // here appropriately - elem_solution(j) = original_solution - numerical_point_h; - *coord = libmesh_real(elem_solution(j)); + context.elem_solution(j) = original_solution - numerical_point_h; + *coord = libmesh_real(context.elem_solution(j)); } - elem_residual.zero(); - ((*time_solver).*(res))(false); + context.elem_residual.zero(); + ((*time_solver).*(res))(false, context); #ifdef DEBUG - libmesh_assert(old_jacobian == elem_jacobian); + libmesh_assert(old_jacobian == context.elem_jacobian); #endif - backwards_residual = elem_residual; + backwards_residual = context.elem_residual; // Take the "plus" side of a central differenced first derivative - elem_solution(j) = original_solution + numerical_jacobian_h; + context.elem_solution(j) = original_solution + numerical_jacobian_h; if (coord) { - elem_solution(j) = original_solution + numerical_point_h; - *coord = libmesh_real(elem_solution(j)); + context.elem_solution(j) = original_solution + numerical_point_h; + *coord = libmesh_real(context.elem_solution(j)); } - elem_residual.zero(); - ((*time_solver).*(res))(false); + context.elem_residual.zero(); + ((*time_solver).*(res))(false, context); #ifdef DEBUG - libmesh_assert(old_jacobian == elem_jacobian); + libmesh_assert(old_jacobian == context.elem_jacobian); #endif - elem_solution(j) = original_solution; + context.elem_solution(j) = original_solution; if (coord) { - *coord = libmesh_real(elem_solution(j)); - for (unsigned int i = 0; i != dof_indices.size(); ++i) + *coord = libmesh_real(context.elem_solution(j)); + for (unsigned int i = 0; i != context.dof_indices.size(); ++i) { numerical_jacobian(i,j) = - (elem_residual(i) - backwards_residual(i)) / + (context.elem_residual(i) - backwards_residual(i)) / 2. / numerical_point_h; } } else { - for (unsigned int i = 0; i != dof_indices.size(); ++i) + for (unsigned int i = 0; i != context.dof_indices.size(); ++i) { numerical_jacobian(i,j) = - (elem_residual(i) - backwards_residual(i)) / + (context.elem_residual(i) - backwards_residual(i)) / 2. / numerical_jacobian_h; } } } - elem_residual = original_residual; - elem_jacobian = numerical_jacobian; + context.elem_residual = original_residual; + context.elem_jacobian = numerical_jacobian; } -void FEMSystem::numerical_elem_jacobian () +void FEMSystem::numerical_elem_jacobian (FEMContext &context) { START_LOG("numerical_elem_jacobian()", "FEMSystem"); - this->numerical_jacobian(&TimeSolver::element_residual); + this->numerical_jacobian(&TimeSolver::element_residual, context); STOP_LOG("numerical_elem_jacobian()", "FEMSystem"); } -void FEMSystem::numerical_side_jacobian () +void FEMSystem::numerical_side_jacobian (FEMContext &context) { START_LOG("numerical_side_jacobian()", "FEMSystem"); - this->numerical_jacobian(&TimeSolver::side_residual); + this->numerical_jacobian(&TimeSolver::side_residual, context); STOP_LOG("numerical_side_jacobian()", "FEMSystem"); } +AutoPtr FEMSystem::build_context () +{ + return AutoPtr(new FEMContext(*this)); +} + + + +void FEMSystem::init_context(DiffContext &c) +{ + Parent::init_context(c); + + FEMContext &context = libmesh_cast_ref(c); + + // Make sure we're prepared to do mass integration + for (unsigned int var = 0; var != this->n_vars(); ++var) + if (_time_evolving[var]) + { + context.element_fe_var[var]->get_JxW(); + context.element_fe_var[var]->get_phi(); + } +} + + + void FEMSystem::time_evolving (unsigned int var) { // Call the parent function Parent::time_evolving(var); - - // Then make sure we're prepared to do mass integration - element_fe_var[var]->get_JxW(); - element_fe_var[var]->get_phi(); } @@ -1264,148 +693,59 @@ const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); + AutoPtr con = this->build_context(); + FEMContext &_femcontext = libmesh_cast_ref(*con); + this->init_context(_femcontext); + // Get the solution's mesh variables from every element for ( ; el != end_el; ++el) { - elem = *el; + _femcontext.elem = *el; // Initialize the per-variable data for elem. for (unsigned int i=0; i != this->n_vars(); ++i) { - dof_map.dof_indices (elem, dof_indices_var[i], i); + dof_map.dof_indices (_femcontext.elem, + _femcontext.dof_indices_var[i], i); } - this->elem_position_get(); - } - - // And make sure the current_local_solution is up to date too - this->update(); -} - - - -void FEMSystem::elem_position_get() -{ - // This is too expensive to call unless we've been asked to move the mesh - libmesh_assert (_mesh_sys != libMesh::invalid_uint); - - // If the coordinate data is in our own system, it's already - // been set up for us - if (_mesh_sys == this->number()) - { - unsigned int n_nodes = elem->n_nodes(); - // For simplicity we demand that mesh coordinates be stored - // in a format that allows a direct copy - libmesh_assert(_mesh_x_var == libMesh::invalid_uint || - (element_fe_var[_mesh_x_var]->get_fe_type().family - == LAGRANGE && - element_fe_var[_mesh_x_var]->get_fe_type().order - == elem->default_order())); - libmesh_assert(_mesh_y_var == libMesh::invalid_uint || - (element_fe_var[_mesh_y_var]->get_fe_type().family - == LAGRANGE && - element_fe_var[_mesh_y_var]->get_fe_type().order - == elem->default_order())); - libmesh_assert(_mesh_z_var == libMesh::invalid_uint || - (element_fe_var[_mesh_z_var]->get_fe_type().family - == LAGRANGE && - element_fe_var[_mesh_z_var]->get_fe_type().order - == elem->default_order())); - - // Get degree of freedom coefficients from point coordinates + _femcontext.elem_position_get(); if (_mesh_x_var != libMesh::invalid_uint) - for (unsigned int i=0; i != n_nodes; ++i) - this->solution->set(dof_indices_var[_mesh_x_var][i], - elem->point(i)(0)); - + this->solution->insert(*_femcontext.elem_subsolutions[_mesh_x_var], + _femcontext.dof_indices_var[_mesh_x_var]); if (_mesh_y_var != libMesh::invalid_uint) - for (unsigned int i=0; i != n_nodes; ++i) - this->solution->set(dof_indices_var[_mesh_y_var][i], - elem->point(i)(1)); - + this->solution->insert(*_femcontext.elem_subsolutions[_mesh_y_var], + _femcontext.dof_indices_var[_mesh_y_var]); if (_mesh_z_var != libMesh::invalid_uint) - for (unsigned int i=0; i != n_nodes; ++i) - this->solution->set(dof_indices_var[_mesh_z_var][i], - elem->point(i)(2)); + this->solution->insert(*_femcontext.elem_subsolutions[_mesh_z_var], + _femcontext.dof_indices_var[_mesh_z_var]); } - // FIXME - If the coordinate data is not in our own system, someone - // had better get around to implementing that... - RHS - else - { - libmesh_not_implemented(); - } -} - - -void FEMSystem::elem_position_set(Real) -{ - // This is too expensive to call unless we've been asked to move the mesh - libmesh_assert (_mesh_sys != libMesh::invalid_uint); - - // If the coordinate data is in our own system, it's already - // been set up for us, and we can ignore our input parameter theta - if (_mesh_sys == this->number()) - { - unsigned int n_nodes = elem->n_nodes(); - // For simplicity we demand that mesh coordinates be stored - // in a format that allows a direct copy - libmesh_assert(_mesh_x_var == libMesh::invalid_uint || - (element_fe_var[_mesh_x_var]->get_fe_type().family - == LAGRANGE && - elem_subsolutions[_mesh_x_var]->size() == n_nodes)); - libmesh_assert(_mesh_y_var == libMesh::invalid_uint || - (element_fe_var[_mesh_y_var]->get_fe_type().family - == LAGRANGE && - elem_subsolutions[_mesh_y_var]->size() == n_nodes)); - libmesh_assert(_mesh_z_var == libMesh::invalid_uint || - (element_fe_var[_mesh_z_var]->get_fe_type().family - == LAGRANGE && - elem_subsolutions[_mesh_z_var]->size() == n_nodes)); - - // Set the new point coordinates - if (_mesh_x_var != libMesh::invalid_uint) - for (unsigned int i=0; i != n_nodes; ++i) - elem->point(i)(0) = - libmesh_real((*elem_subsolutions[_mesh_x_var])(i)); - - if (_mesh_y_var != libMesh::invalid_uint) - for (unsigned int i=0; i != n_nodes; ++i) - elem->point(i)(1) = - libmesh_real((*elem_subsolutions[_mesh_y_var])(i)); - - if (_mesh_z_var != libMesh::invalid_uint) - for (unsigned int i=0; i != n_nodes; ++i) - elem->point(i)(2) = - libmesh_real((*elem_subsolutions[_mesh_z_var])(i)); - } - // FIXME - If the coordinate data is not in our own system, someone - // had better get around to implementing that... - RHS - else - { - libmesh_not_implemented(); - } + // And make sure the current_local_solution is up to date too + this->update(); } - -bool FEMSystem::eulerian_residual (bool request_jacobian) +bool FEMSystem::eulerian_residual (bool request_jacobian, + DiffContext &c) { // Only calculate a mesh movement residual if it's necessary if (_mesh_sys == libMesh::invalid_uint) return request_jacobian; + FEMContext &context = libmesh_cast_ref(c); + // This function only supports fully coupled mesh motion for now libmesh_assert(_mesh_sys == this->number()); - unsigned int n_qpoints = element_qrule->n_points(); + unsigned int n_qpoints = context.element_qrule->n_points(); const unsigned int n_x_dofs = (_mesh_x_var == libMesh::invalid_uint) ? - 0 : dof_indices_var[_mesh_x_var].size(); + 0 : context.dof_indices_var[_mesh_x_var].size(); const unsigned int n_y_dofs = (_mesh_y_var == libMesh::invalid_uint) ? - 0 : dof_indices_var[_mesh_y_var].size(); + 0 : context.dof_indices_var[_mesh_y_var].size(); const unsigned int n_z_dofs = (_mesh_z_var == libMesh::invalid_uint) ? - 0 : dof_indices_var[_mesh_z_var].size(); + 0 : context.dof_indices_var[_mesh_z_var].size(); const unsigned int mesh_xyz_var = n_x_dofs ? _mesh_x_var : (n_y_dofs ? _mesh_y_var : @@ -1416,15 +756,15 @@ // at least one coordinate, and we'd better have the same // FE type for all coordinates we are in charge of libmesh_assert(mesh_xyz_var != libMesh::invalid_uint); - libmesh_assert(!n_x_dofs || element_fe_var[_mesh_x_var] == - element_fe_var[mesh_xyz_var]); - libmesh_assert(!n_y_dofs || element_fe_var[_mesh_y_var] == - element_fe_var[mesh_xyz_var]); - libmesh_assert(!n_z_dofs || element_fe_var[_mesh_z_var] == - element_fe_var[mesh_xyz_var]); + libmesh_assert(!n_x_dofs || context.element_fe_var[_mesh_x_var] == + context.element_fe_var[mesh_xyz_var]); + libmesh_assert(!n_y_dofs || context.element_fe_var[_mesh_y_var] == + context.element_fe_var[mesh_xyz_var]); + libmesh_assert(!n_z_dofs || context.element_fe_var[_mesh_z_var] == + context.element_fe_var[mesh_xyz_var]); const std::vector > &psi = - element_fe_var[mesh_xyz_var]->get_phi(); + context.element_fe_var[mesh_xyz_var]->get_phi(); for (unsigned int var = 0; var != this->n_vars(); ++var) { @@ -1453,25 +793,25 @@ return request_jacobian; const std::vector &JxW = - element_fe_var[var]->get_JxW(); + context.element_fe_var[var]->get_JxW(); const std::vector > &phi = - element_fe_var[var]->get_phi(); + context.element_fe_var[var]->get_phi(); const std::vector > &dphi = - element_fe_var[var]->get_dphi(); + context.element_fe_var[var]->get_dphi(); - const unsigned int n_u_dofs = dof_indices_var[var].size(); + const unsigned int n_u_dofs = context.dof_indices_var[var].size(); - DenseSubVector &Fu = *elem_subresiduals[var]; - DenseSubMatrix &Kuu = *elem_subjacobians[var][var]; + DenseSubVector &Fu = *context.elem_subresiduals[var]; + DenseSubMatrix &Kuu = *context.elem_subjacobians[var][var]; DenseSubMatrix *Kux = n_x_dofs ? - elem_subjacobians[var][_mesh_x_var] : NULL; + context.elem_subjacobians[var][_mesh_x_var] : NULL; DenseSubMatrix *Kuy = n_y_dofs ? - elem_subjacobians[var][_mesh_y_var] : NULL; + context.elem_subjacobians[var][_mesh_y_var] : NULL; DenseSubMatrix *Kuz = n_z_dofs ? - elem_subjacobians[var][_mesh_z_var] : NULL; + context.elem_subjacobians[var][_mesh_z_var] : NULL; std::vector delta_x(n_x_dofs, 0.); std::vector delta_y(n_y_dofs, 0.); @@ -1479,28 +819,28 @@ for (unsigned int i = 0; i != n_x_dofs; ++i) { - unsigned int j = dof_indices_var[_mesh_x_var][i]; + unsigned int j = context.dof_indices_var[_mesh_x_var][i]; delta_x[i] = libmesh_real(this->current_solution(j)) - libmesh_real(unsteady->old_nonlinear_solution(j)); } for (unsigned int i = 0; i != n_y_dofs; ++i) { - unsigned int j = dof_indices_var[_mesh_y_var][i]; + unsigned int j = context.dof_indices_var[_mesh_y_var][i]; delta_y[i] = libmesh_real(this->current_solution(j)) - libmesh_real(unsteady->old_nonlinear_solution(j)); } for (unsigned int i = 0; i != n_z_dofs; ++i) { - unsigned int j = dof_indices_var[_mesh_z_var][i]; + unsigned int j = context.dof_indices_var[_mesh_z_var][i]; delta_z[i] = libmesh_real(this->current_solution(j)) - libmesh_real(unsteady->old_nonlinear_solution(j)); } for (unsigned int qp = 0; qp != n_qpoints; ++qp) { - Gradient grad_u = interior_gradient(var, qp); + Gradient grad_u = context.interior_gradient(var, qp); RealGradient convection(0.); for (unsigned int i = 0; i != n_x_dofs; ++i) @@ -1543,36 +883,39 @@ -bool FEMSystem::mass_residual (bool request_jacobian) +bool FEMSystem::mass_residual (bool request_jacobian, + DiffContext &c) { - unsigned int n_qpoints = element_qrule->n_points(); + FEMContext &context = libmesh_cast_ref(c); + unsigned int n_qpoints = context.element_qrule->n_points(); + for (unsigned int var = 0; var != this->n_vars(); ++var) { if (!_time_evolving[var]) continue; const std::vector &JxW = - element_fe_var[var]->get_JxW(); + context.element_fe_var[var]->get_JxW(); const std::vector > &phi = - element_fe_var[var]->get_phi(); + context.element_fe_var[var]->get_phi(); - const unsigned int n_dofs = dof_indices_var[var].size(); + const unsigned int n_dofs = context.dof_indices_var[var].size(); - DenseSubVector &Fu = *elem_subresiduals[var]; - DenseSubMatrix &Kuu = *elem_subjacobians[var][var]; + DenseSubVector &Fu = *context.elem_subresiduals[var]; + DenseSubMatrix &Kuu = *context.elem_subjacobians[var][var]; for (unsigned int qp = 0; qp != n_qpoints; ++qp) { - Number u = interior_value(var, qp); + Number u = context.interior_value(var, qp); Number JxWxU = JxW[qp] * u; for (unsigned int i = 0; i != n_dofs; ++i) { Fu(i) += JxWxU * phi[i][qp]; - if (request_jacobian && elem_solution_derivative) + if (request_jacobian && context.elem_solution_derivative) { - libmesh_assert (elem_solution_derivative == 1.0); + libmesh_assert (context.elem_solution_derivative == 1.0); Number JxWxPhiI = JxW[qp] * phi[i][qp]; Kuu(i,i) += JxWxPhiI * phi[i][qp]; Index: src/solvers/diff_context.C =================================================================== --- src/solvers/diff_context.C (revision 0) +++ src/solvers/diff_context.C (revision 0) @@ -0,0 +1,74 @@ +// $Id: diff_system.C 3118 2008-10-27 15:41:13Z jwpeterson $ + +// The libMesh Finite Element Library. +// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// 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 02111-1307 USA + + +#include "diff_context.h" +#include "diff_system.h" + + + +DiffContext::DiffContext (const DifferentiableSystem& sys) : + time(sys.time), + elem_solution_derivative(1.), + fixed_solution_derivative(0.), + dof_indices_var(sys.n_vars()) +{ + // Finally initialize solution/residual/jacobian data structures + unsigned int n_vars = sys.n_vars(); + + elem_subsolutions.reserve(n_vars); + elem_subresiduals.reserve(n_vars); + elem_subjacobians.resize(n_vars); + if (sys.use_fixed_solution) + elem_fixed_subsolutions.reserve(n_vars); + + for (unsigned int i=0; i != n_vars; ++i) + { + elem_subsolutions.push_back(new DenseSubVector(elem_solution)); + elem_subresiduals.push_back(new DenseSubVector(elem_residual)); + elem_subjacobians[i].reserve(n_vars); + + if (sys.use_fixed_solution) + elem_fixed_subsolutions.push_back + (new DenseSubVector(elem_fixed_solution)); + + for (unsigned int j=0; j != n_vars; ++j) + elem_subjacobians[i].push_back + (new DenseSubMatrix(elem_jacobian)); + } +} + + + +DiffContext::~DiffContext () +{ + for (unsigned int i=0; i != elem_subsolutions.size(); ++i) + { + delete elem_subsolutions[i]; + delete elem_subresiduals[i]; + if (!elem_fixed_subsolutions.empty()) + delete elem_fixed_subsolutions[i]; + + for (unsigned int j=0; j != elem_subjacobians[i].size(); ++j) + delete elem_subjacobians[i][j]; + } +} + + + Index: src/solvers/steady_solver.C =================================================================== --- src/solvers/steady_solver.C (revision 3205) +++ src/solvers/steady_solver.C (working copy) @@ -30,55 +30,57 @@ -bool SteadySolver::element_residual(bool request_jacobian) +bool SteadySolver::element_residual(bool request_jacobian, + DiffContext& context) { // If a fixed solution is requested, it will just be the current // solution if (_system.use_fixed_solution) { - _system.elem_fixed_solution = _system.elem_solution; - _system.fixed_solution_derivative = 1.0; + context.elem_fixed_solution = context.elem_solution; + context.fixed_solution_derivative = 1.0; } bool jacobian_computed = - _system.element_time_derivative(request_jacobian); + _system.element_time_derivative(request_jacobian, context); // The user shouldn't compute a jacobian unless requested libmesh_assert(request_jacobian || !jacobian_computed); bool jacobian_computed2 = - _system.element_constraint(jacobian_computed); + _system.element_constraint(jacobian_computed, context); // The user shouldn't compute a jacobian unless requested libmesh_assert (jacobian_computed || !jacobian_computed2); - return jacobian_computed && jacobian_computed2; + return jacobian_computed2; } -bool SteadySolver::side_residual(bool request_jacobian) +bool SteadySolver::side_residual(bool request_jacobian, + DiffContext& context) { // If a fixed solution is requested, it will just be the current // solution if (_system.use_fixed_solution) { - _system.elem_fixed_solution = _system.elem_solution; - _system.fixed_solution_derivative = 1.0; + context.elem_fixed_solution = context.elem_solution; + context.fixed_solution_derivative = 1.0; } bool jacobian_computed = - _system.side_time_derivative(request_jacobian); + _system.side_time_derivative(request_jacobian, context); // The user shouldn't compute a jacobian unless requested libmesh_assert (request_jacobian || !jacobian_computed); bool jacobian_computed2 = - _system.side_constraint(jacobian_computed); + _system.side_constraint(jacobian_computed, context); // The user shouldn't compute a jacobian unless requested libmesh_assert (jacobian_computed || !jacobian_computed2); - return jacobian_computed && jacobian_computed2; + return jacobian_computed2; } Index: src/solvers/euler2_solver.C =================================================================== --- src/solvers/euler2_solver.C (revision 3205) +++ src/solvers/euler2_solver.C (working copy) @@ -47,64 +47,65 @@ -bool Euler2Solver::element_residual (bool request_jacobian) +bool Euler2Solver::element_residual (bool request_jacobian, + DiffContext &context) { - unsigned int n_dofs = _system.elem_solution.size(); + unsigned int n_dofs = context.elem_solution.size(); // Local nonlinear solution at old timestep DenseVector old_elem_solution(n_dofs); for (unsigned int i=0; i != n_dofs; ++i) old_elem_solution(i) = - old_nonlinear_solution(_system.dof_indices[i]); + old_nonlinear_solution(context.dof_indices[i]); // We evaluate mass_residual with the change in solution // to get the mass matrix, reusing old_elem_solution to hold that // delta_solution. - DenseVector delta_elem_solution(_system.elem_solution); + DenseVector delta_elem_solution(context.elem_solution); delta_elem_solution -= old_elem_solution; // Our first evaluations are at the true elem_solution - _system.elem_solution_derivative = 1.0; + context.elem_solution_derivative = 1.0; // If a fixed solution is requested, we'll use the elem_solution // at the new timestep if (_system.use_fixed_solution) - _system.elem_fixed_solution = _system.elem_solution; + context.elem_fixed_solution = context.elem_solution; - _system.fixed_solution_derivative = 1.0; + context.fixed_solution_derivative = 1.0; // We're going to compute just the change in elem_residual // (and possibly elem_jacobian), then add back the old values - DenseVector total_elem_residual(_system.elem_residual); + DenseVector total_elem_residual(context.elem_residual); DenseMatrix old_elem_jacobian, total_elem_jacobian; if (request_jacobian) { - old_elem_jacobian = _system.elem_jacobian; - total_elem_jacobian = _system.elem_jacobian; - _system.elem_jacobian.zero(); + old_elem_jacobian = context.elem_jacobian; + total_elem_jacobian = context.elem_jacobian; + context.elem_jacobian.zero(); } - _system.elem_residual.zero(); + context.elem_residual.zero(); // First, evaluate time derivative at the new timestep. // The element should already be in the proper place // even for a moving mesh problem. bool jacobian_computed = - _system.element_time_derivative(request_jacobian); + _system.element_time_derivative(request_jacobian, context); // For a moving mesh problem we may need the pseudoconvection term too jacobian_computed = - _system.eulerian_residual(jacobian_computed) && jacobian_computed; + _system.eulerian_residual(jacobian_computed, context) && jacobian_computed; // Scale the new time-dependent residual and jacobian correctly - _system.elem_residual *= (theta * _system.deltat); - total_elem_residual += _system.elem_residual; - _system.elem_residual.zero(); + context.elem_residual *= (theta * _system.deltat); + total_elem_residual += context.elem_residual; + context.elem_residual.zero(); if (jacobian_computed) { - _system.elem_jacobian *= (theta * _system.deltat); - total_elem_jacobian += _system.elem_jacobian; - _system.elem_jacobian.zero(); + context.elem_jacobian *= (theta * _system.deltat); + total_elem_jacobian += context.elem_jacobian; + context.elem_jacobian.zero(); } // Next, evaluate the mass residual at the new timestep, @@ -115,102 +116,104 @@ // mass_residual functions // Move elem_->delta_, delta_->elem_ - _system.elem_solution.swap(delta_elem_solution); + context.elem_solution.swap(delta_elem_solution); - jacobian_computed = _system.mass_residual(jacobian_computed) && + jacobian_computed = _system.mass_residual(jacobian_computed, context) && jacobian_computed; - _system.elem_residual *= -theta; - total_elem_residual += _system.elem_residual; - _system.elem_residual.zero(); + context.elem_residual *= -theta; + total_elem_residual += context.elem_residual; + context.elem_residual.zero(); if (jacobian_computed) { // The minus sign trick here is to avoid using a non-1 // elem_solution_derivative: - _system.elem_jacobian *= -theta; - total_elem_jacobian += _system.elem_jacobian; - _system.elem_jacobian.zero(); + context.elem_jacobian *= -theta; + total_elem_jacobian += context.elem_jacobian; + context.elem_jacobian.zero(); } // Add the time-dependent term for the old solution // Make sure elem_solution is set up for elem_reinit to use // Move delta_->old_, old_->elem_ - _system.elem_solution.swap(old_elem_solution); + context.elem_solution.swap(old_elem_solution); // Move the mesh into place first if necessary - _system.elem_reinit(0.); + context.elem_reinit(0.); if (_system.use_fixed_solution) { - _system.elem_solution_derivative = 0.0; - jacobian_computed =_system.element_time_derivative(jacobian_computed) && + context.elem_solution_derivative = 0.0; + jacobian_computed = + _system.element_time_derivative(jacobian_computed, context) && jacobian_computed; - jacobian_computed = _system.eulerian_residual(jacobian_computed) && + jacobian_computed = + _system.eulerian_residual(jacobian_computed, context) && jacobian_computed; - _system.elem_solution_derivative = 1.0; - _system.elem_residual *= ((1. - theta) * _system.deltat); - total_elem_residual += _system.elem_residual; + context.elem_solution_derivative = 1.0; + context.elem_residual *= ((1. - theta) * _system.deltat); + total_elem_residual += context.elem_residual; if (jacobian_computed) { - _system.elem_jacobian *= ((1. - theta) * _system.deltat); - total_elem_jacobian += _system.elem_jacobian; - _system.elem_jacobian.zero(); + context.elem_jacobian *= ((1. - theta) * _system.deltat); + total_elem_jacobian += context.elem_jacobian; + context.elem_jacobian.zero(); } } else { // FIXME - we should detect if element_time_derivative() edits // elem_jacobian and lies about it! - _system.element_time_derivative(false); - _system.eulerian_residual(false); - _system.elem_residual *= ((1. - theta) * _system.deltat); - total_elem_residual += _system.elem_residual; + _system.element_time_derivative(false, context); + _system.eulerian_residual(false, context); + context.elem_residual *= ((1. - theta) * _system.deltat); + total_elem_residual += context.elem_residual; } - _system.elem_residual.zero(); + context.elem_residual.zero(); // Add the mass residual term for the old solution // Move old_->old_, delta_->elem_ - _system.elem_solution.swap(old_elem_solution); + context.elem_solution.swap(old_elem_solution); - jacobian_computed = _system.mass_residual(jacobian_computed) && + jacobian_computed = _system.mass_residual(jacobian_computed, context) && jacobian_computed; - _system.elem_residual *= -(1. - theta); - total_elem_residual += _system.elem_residual; - _system.elem_residual.zero(); + context.elem_residual *= -(1. - theta); + total_elem_residual += context.elem_residual; + context.elem_residual.zero(); if (jacobian_computed) { // The minus sign trick here is to avoid using a non-1 // *_solution_derivative: - _system.elem_jacobian *= -(1. - theta); - total_elem_jacobian += _system.elem_jacobian; - _system.elem_jacobian.zero(); + context.elem_jacobian *= -(1. - theta); + total_elem_jacobian += context.elem_jacobian; + context.elem_jacobian.zero(); } // Restore the elem_solution // Move elem_->elem_, delta_->delta_ - _system.elem_solution.swap(delta_elem_solution); + context.elem_solution.swap(delta_elem_solution); // Restore the elem position if necessary - _system.elem_reinit(1.); + context.elem_reinit(1.); // Add the constraint term - jacobian_computed = _system.element_constraint(jacobian_computed) && + jacobian_computed = _system.element_constraint(jacobian_computed, context) && jacobian_computed; // Add back the previously accumulated residual and jacobian - _system.elem_residual += total_elem_residual; + context.elem_residual += total_elem_residual; if (request_jacobian) { if (jacobian_computed) - _system.elem_jacobian += total_elem_jacobian; + context.elem_jacobian += total_elem_jacobian; else - _system.elem_jacobian.swap(old_elem_jacobian); + context.elem_jacobian.swap(old_elem_jacobian); } return jacobian_computed; @@ -218,60 +221,61 @@ -bool Euler2Solver::side_residual (bool request_jacobian) +bool Euler2Solver::side_residual (bool request_jacobian, + DiffContext &context) { - unsigned int n_dofs = _system.elem_solution.size(); + unsigned int n_dofs = context.elem_solution.size(); // Local nonlinear solution at old timestep DenseVector old_elem_solution(n_dofs); for (unsigned int i=0; i != n_dofs; ++i) old_elem_solution(i) = - old_nonlinear_solution(_system.dof_indices[i]); + old_nonlinear_solution(context.dof_indices[i]); // We evaluate mass_residual with the change in solution // to get the mass matrix, reusing old_elem_solution to hold that // delta_solution. - DenseVector delta_elem_solution(_system.elem_solution); + DenseVector delta_elem_solution(context.elem_solution); delta_elem_solution -= old_elem_solution; // Our first evaluations are at the true elem_solution - _system.elem_solution_derivative = 1.0; + context.elem_solution_derivative = 1.0; // If a fixed solution is requested, we'll use the elem_solution // at the new timestep if (_system.use_fixed_solution) - _system.elem_fixed_solution = _system.elem_solution; + context.elem_fixed_solution = context.elem_solution; - _system.fixed_solution_derivative = 1.0; + context.fixed_solution_derivative = 1.0; // We're going to compute just the change in elem_residual // (and possibly elem_jacobian), then add back the old values - DenseVector total_elem_residual(_system.elem_residual); + DenseVector total_elem_residual(context.elem_residual); DenseMatrix old_elem_jacobian, total_elem_jacobian; if (request_jacobian) { - old_elem_jacobian = _system.elem_jacobian; - total_elem_jacobian = _system.elem_jacobian; - _system.elem_jacobian.zero(); + old_elem_jacobian = context.elem_jacobian; + total_elem_jacobian = context.elem_jacobian; + context.elem_jacobian.zero(); } - _system.elem_residual.zero(); + context.elem_residual.zero(); // First, evaluate time derivative at the new timestep. // The element should already be in the proper place // even for a moving mesh problem. bool jacobian_computed = - _system.side_time_derivative(request_jacobian); + _system.side_time_derivative(request_jacobian, context); // Scale the time-dependent residual and jacobian correctly - _system.elem_residual *= (theta * _system.deltat); - total_elem_residual += _system.elem_residual; - _system.elem_residual.zero(); + context.elem_residual *= (theta * _system.deltat); + total_elem_residual += context.elem_residual; + context.elem_residual.zero(); if (jacobian_computed) { - _system.elem_jacobian *= (theta * _system.deltat); - total_elem_jacobian += _system.elem_jacobian; - _system.elem_jacobian.zero(); + context.elem_jacobian *= (theta * _system.deltat); + total_elem_jacobian += context.elem_jacobian; + context.elem_jacobian.zero(); } // Next, evaluate the mass residual at the new timestep, @@ -282,99 +286,101 @@ // mass_residual functions // Move elem_->delta_, delta_->elem_ - _system.elem_solution.swap(delta_elem_solution); + context.elem_solution.swap(delta_elem_solution); - jacobian_computed = _system.side_mass_residual(jacobian_computed) && + jacobian_computed = _system.side_mass_residual(jacobian_computed, context) && jacobian_computed; - _system.elem_residual *= -theta; - total_elem_residual += _system.elem_residual; - _system.elem_residual.zero(); + context.elem_residual *= -theta; + total_elem_residual += context.elem_residual; + context.elem_residual.zero(); if (jacobian_computed) { // The minus sign trick here is to avoid using a non-1 // elem_solution_derivative: - _system.elem_jacobian *= -theta; - total_elem_jacobian += _system.elem_jacobian; - _system.elem_jacobian.zero(); + context.elem_jacobian *= -theta; + total_elem_jacobian += context.elem_jacobian; + context.elem_jacobian.zero(); } // Add the time-dependent term for the old solution // Make sure elem_solution is set up for elem_reinit to use // Move delta_->old_, old_->elem_ - _system.elem_solution.swap(old_elem_solution); + context.elem_solution.swap(old_elem_solution); // Move the mesh into place first if necessary - _system.elem_side_reinit(0.); + context.elem_side_reinit(0.); if (_system.use_fixed_solution) { - _system.elem_solution_derivative = 0.0; - jacobian_computed = _system.side_time_derivative(jacobian_computed) && - jacobian_computed; - _system.elem_solution_derivative = 1.0; - _system.elem_residual *= ((1. - theta) * _system.deltat); - total_elem_residual += _system.elem_residual; + context.elem_solution_derivative = 0.0; + jacobian_computed = + _system.side_time_derivative(jacobian_computed, context) && + jacobian_computed; + context.elem_solution_derivative = 1.0; + context.elem_residual *= ((1. - theta) * _system.deltat); + total_elem_residual += context.elem_residual; if (jacobian_computed) { - _system.elem_jacobian *= ((1. - theta) * _system.deltat); - total_elem_jacobian += _system.elem_jacobian; - _system.elem_jacobian.zero(); + context.elem_jacobian *= ((1. - theta) * _system.deltat); + total_elem_jacobian += context.elem_jacobian; + context.elem_jacobian.zero(); } } else { // FIXME - we should detect if side_time_derivative() edits // elem_jacobian and lies about it! - _system.side_time_derivative(false); - _system.elem_residual *= ((1. - theta) * _system.deltat); - total_elem_residual += _system.elem_residual; + _system.side_time_derivative(false, context); + context.elem_residual *= ((1. - theta) * _system.deltat); + total_elem_residual += context.elem_residual; } - _system.elem_residual.zero(); + context.elem_residual.zero(); // Add the mass residual term for the old solution // Move old_->old_, delta_->elem_ - _system.elem_solution.swap(old_elem_solution); + context.elem_solution.swap(old_elem_solution); - jacobian_computed = _system.side_mass_residual(jacobian_computed) && - jacobian_computed; + jacobian_computed = + _system.side_mass_residual(jacobian_computed, context) && + jacobian_computed; - _system.elem_residual *= -(1. - theta); - total_elem_residual += _system.elem_residual; - _system.elem_residual.zero(); + context.elem_residual *= -(1. - theta); + total_elem_residual += context.elem_residual; + context.elem_residual.zero(); if (jacobian_computed) { // The minus sign trick here is to avoid using a non-1 // *_solution_derivative: - _system.elem_jacobian *= -(1. - theta); - total_elem_jacobian += _system.elem_jacobian; - _system.elem_jacobian.zero(); + context.elem_jacobian *= -(1. - theta); + total_elem_jacobian += context.elem_jacobian; + context.elem_jacobian.zero(); } // Restore the elem_solution // Move elem_->elem_, delta_->delta_ - _system.elem_solution.swap(delta_elem_solution); + context.elem_solution.swap(delta_elem_solution); // Restore the elem position if necessary - _system.elem_side_reinit(1.); + context.elem_side_reinit(1.); // Add the constraint term - jacobian_computed = _system.side_constraint(jacobian_computed) && + jacobian_computed = _system.side_constraint(jacobian_computed, context) && jacobian_computed; // Add back the previously accumulated residual and jacobian - _system.elem_residual += total_elem_residual; + context.elem_residual += total_elem_residual; if (request_jacobian) { if (jacobian_computed) - _system.elem_jacobian += total_elem_jacobian; + context.elem_jacobian += total_elem_jacobian; else - _system.elem_jacobian.swap(old_elem_jacobian); + context.elem_jacobian.swap(old_elem_jacobian); } return jacobian_computed; Index: .depend =================================================================== --- .depend (revision 3205) +++ .depend (working copy) @@ -45,6 +45,7 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h\ @@ -98,6 +99,7 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/solvers/system.h\ include/solvers/system_norm.h\ include/utils/compare_types.h\ @@ -209,6 +211,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -279,6 +282,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ @@ -343,6 +347,7 @@ include/numerics/type_vector.h\ include/numerics/vector_value.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/system.h\ include/solvers/system_norm.h\ include/utils/compare_types.h\ @@ -401,6 +406,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/quadrature/quadrature_gauss.h\ include/solvers/system.h\ @@ -457,6 +463,7 @@ include/numerics/type_vector.h\ include/numerics/vector_value.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/quadrature/quadrature_gauss.h\ include/solvers/system.h\ @@ -519,6 +526,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/quadrature/quadrature_gauss.h\ include/quadrature/quadrature_grid.h\ @@ -581,6 +589,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ @@ -688,6 +697,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/quadrature/quadrature_gauss.h\ include/utils/compare_types.h\ @@ -1069,6 +1079,7 @@ include/mesh/mesh_base.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -1429,6 +1440,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -1503,6 +1515,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -3552,6 +3565,7 @@ include/numerics/type_vector.h\ include/numerics/vector_value.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/quadrature/quadrature_gauss.h\ include/utils/compare_types.h\ @@ -4109,6 +4123,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/mapvector.h\ include/utils/o_string_stream.h\ @@ -4145,6 +4160,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/mapvector.h\ include/utils/o_string_stream.h\ @@ -4185,6 +4201,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/mapvector.h\ include/utils/o_string_stream.h\ @@ -4229,6 +4246,7 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/system.h\ include/solvers/system_norm.h\ include/utils/compare_types.h\ @@ -4272,6 +4290,7 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/system.h\ include/solvers/system_norm.h\ include/utils/compare_types.h\ @@ -4307,6 +4326,7 @@ include/mesh/mesh_output.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -4341,6 +4361,7 @@ include/mesh/mesh_output.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -4383,6 +4404,7 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -4426,6 +4448,7 @@ include/mesh/mesh_output.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -4470,6 +4493,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/mapvector.h\ include/utils/o_string_stream.h\ @@ -4517,6 +4541,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_f_stream.h\ include/utils/o_string_stream.h\ @@ -4554,6 +4579,7 @@ include/mesh/mesh_input.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -4586,6 +4612,7 @@ include/mesh/mesh_output.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -4659,6 +4686,7 @@ include/parallel/parallel.h\ include/parallel/parallel_ghost_sync.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/location_maps.h\ include/utils/mapvector.h\ @@ -4699,6 +4727,7 @@ include/parallel/parallel.h\ include/parallel/parallel_sort.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -4731,6 +4760,7 @@ include/mesh/mesh_data.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -4790,6 +4820,7 @@ include/mesh/mesh_data.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h\ @@ -4854,6 +4885,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/location_maps.h\ include/utils/o_string_stream.h\ @@ -4895,6 +4927,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/location_maps.h\ include/utils/o_string_stream.h\ @@ -4936,6 +4969,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -4986,6 +5020,7 @@ include/parallel/parallel.h\ include/parallel/parallel_ghost_sync.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/error_vector.h\ include/utils/location_maps.h\ @@ -5022,6 +5057,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/error_vector.h\ include/utils/location_maps.h\ @@ -5058,6 +5094,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/location_maps.h\ include/utils/o_string_stream.h\ @@ -5095,6 +5132,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5129,6 +5167,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h\ @@ -5168,6 +5207,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5210,6 +5250,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/location_maps.h\ include/utils/mapvector.h\ @@ -5253,6 +5294,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5292,6 +5334,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/mapvector.h\ include/utils/o_string_stream.h\ @@ -5319,6 +5362,7 @@ include/mesh/mesh_base.h\ include/mesh/nemesis_io_helper.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h src/mesh/off_io.$(obj-suffix):\ @@ -5353,6 +5397,7 @@ include/mesh/off_io.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5424,6 +5469,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/mapvector.h\ include/utils/o_string_stream.h\ @@ -5488,6 +5534,7 @@ include/mesh/postscript_io.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5555,6 +5602,7 @@ include/mesh/unstructured_mesh.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5587,6 +5635,7 @@ include/mesh/tecplot_io.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5626,6 +5675,7 @@ include/mesh/tetgen_io.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5670,6 +5720,7 @@ include/mesh/ucd_io.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5725,6 +5776,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5776,6 +5828,7 @@ include/mesh/unv_io.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5832,6 +5885,7 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -5872,6 +5926,7 @@ include/mesh/xdr_head.h\ include/mesh/xdr_mgf.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/o_f_stream.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5942,6 +5997,7 @@ include/mesh/xdr_mgf.h\ include/mesh/xdr_mhead.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/o_f_stream.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5969,6 +6025,7 @@ include/mesh/mesh_output.h\ include/mesh/xdr_mgf.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/o_f_stream.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -5999,6 +6056,7 @@ include/mesh/xdr_shead.h\ include/mesh/xdr_soln.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/o_f_stream.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -6127,6 +6185,7 @@ include/numerics/type_vector.h\ include/numerics/vector_value.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/system.h\ include/solvers/system_norm.h\ include/utils/compare_types.h\ @@ -6147,6 +6206,7 @@ include/base/reference_counted_object.h\ include/base/reference_counter.h\ include/enums/enum_solver_package.h\ + include/numerics/dense_subvector.h\ include/numerics/dense_vector.h\ include/numerics/dense_vector_base.h\ include/numerics/distributed_vector.h\ @@ -6242,6 +6302,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/solvers/system.h\ include/solvers/system_norm.h\ @@ -6308,6 +6369,7 @@ include/numerics/type_vector.h\ include/numerics/vector_value.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/solvers/system.h\ include/solvers/system_norm.h\ @@ -6332,8 +6394,6 @@ include/enums/enum_solver_type.h\ include/numerics/dense_matrix.h\ include/numerics/dense_matrix_base.h\ - include/numerics/dense_vector.h\ - include/numerics/dense_vector_base.h\ include/numerics/laspack_linear_solver.h\ include/numerics/laspack_matrix.h\ include/numerics/laspack_vector.h\ @@ -6380,6 +6440,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -6398,6 +6459,7 @@ include/enums/enum_solver_package.h\ include/numerics/dense_matrix.h\ include/numerics/dense_matrix_base.h\ + include/numerics/dense_subvector.h\ include/numerics/dense_vector.h\ include/numerics/dense_vector_base.h\ include/numerics/laspack_matrix.h\ @@ -6426,8 +6488,6 @@ include/enums/enum_solver_type.h\ include/numerics/dense_matrix.h\ include/numerics/dense_matrix_base.h\ - include/numerics/dense_vector.h\ - include/numerics/dense_vector_base.h\ include/numerics/laspack_linear_solver.h\ include/numerics/laspack_matrix.h\ include/numerics/laspack_vector.h\ @@ -6496,6 +6556,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -6547,8 +6608,6 @@ include/base/reference_counted_object.h\ include/base/reference_counter.h\ include/enums/enum_solver_package.h\ - include/numerics/dense_vector.h\ - include/numerics/dense_vector_base.h\ include/numerics/distributed_vector.h\ include/numerics/laspack_vector.h\ include/numerics/numeric_vector.h\ @@ -6558,7 +6617,6 @@ include/numerics/trilinos_epetra_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ - include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h src/numerics/petsc_linear_solver.$(obj-suffix):\ @@ -6628,6 +6686,7 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -6676,6 +6735,7 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/solvers/explicit_system.h\ include/solvers/implicit_system.h\ include/solvers/nonlinear_implicit_system.h\ @@ -6697,6 +6757,7 @@ include/base/reference_counted_object.h\ include/base/reference_counter.h\ include/enums/enum_solver_package.h\ + include/numerics/dense_subvector.h\ include/numerics/dense_vector.h\ include/numerics/dense_vector_base.h\ include/numerics/numeric_vector.h\ @@ -6776,6 +6837,7 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -6899,6 +6961,7 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h @@ -6915,12 +6978,12 @@ include/base/reference_counted_object.h\ include/base/reference_counter.h\ include/enums/enum_solver_package.h\ + include/numerics/dense_subvector.h\ include/numerics/dense_vector.h\ include/numerics/dense_vector_base.h\ include/numerics/numeric_vector.h\ include/numerics/trilinos_epetra_vector.h\ include/parallel/parallel.h\ - include/parallel/parallel.h\ include/parallel/threads.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ @@ -6964,6 +7027,7 @@ include/numerics/trilinos_nox_nonlinear_solver.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/explicit_system.h\ include/solvers/implicit_system.h\ include/solvers/nonlinear_implicit_system.h\ @@ -7996,7 +8060,9 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/adaptive_time_solver.h\ + include/solvers/diff_context.h\ include/solvers/diff_system.h\ include/solvers/euler_solver.h\ include/solvers/explicit_system.h\ @@ -8057,10 +8123,13 @@ include/numerics/type_tensor.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/continuation_system.h\ + include/solvers/diff_context.h\ include/solvers/diff_solver.h\ include/solvers/diff_system.h\ include/solvers/explicit_system.h\ + include/solvers/fem_context.h\ include/solvers/fem_system.h\ include/solvers/implicit_system.h\ include/solvers/linear_implicit_system.h\ @@ -8072,6 +8141,58 @@ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h +src/solvers/diff_context.$(obj-suffix):\ + src/solvers/diff_context.C\ + include/base/auto_ptr.h\ + include/base/dof_object.h\ + include/base/libmesh.h\ + include/base/libmesh_base.h\ + include/base/libmesh_common.h\ + include/base/libmesh_config.h\ + include/base/libmesh_exceptions.h\ + include/base/libmesh_logging.h\ + include/base/multi_predicates.h\ + include/base/print_trace.h\ + include/base/reference_counted_object.h\ + include/base/reference_counter.h\ + include/base/single_predicates.h\ + include/base/variant_filter_iterator.h\ + include/enums/enum_elem_quality.h\ + include/enums/enum_elem_type.h\ + include/enums/enum_fe_family.h\ + include/enums/enum_inf_map_type.h\ + include/enums/enum_io_package.h\ + include/enums/enum_norm_type.h\ + include/enums/enum_order.h\ + include/enums/enum_solver_package.h\ + include/enums/enum_xdr_mode.h\ + include/fe/fe_type.h\ + include/geom/elem.h\ + include/geom/elem_range.h\ + include/geom/node.h\ + include/geom/point.h\ + include/geom/stored_range.h\ + include/mesh/mesh_base.h\ + include/numerics/dense_matrix.h\ + include/numerics/dense_matrix_base.h\ + include/numerics/dense_submatrix.h\ + include/numerics/dense_subvector.h\ + include/numerics/dense_vector.h\ + include/numerics/dense_vector_base.h\ + include/numerics/type_vector.h\ + include/parallel/threads.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ + include/solvers/diff_system.h\ + include/solvers/explicit_system.h\ + include/solvers/implicit_system.h\ + include/solvers/linear_implicit_system.h\ + include/solvers/system.h\ + include/solvers/system_norm.h\ + include/solvers/transient_system.h\ + include/utils/compare_types.h\ + include/utils/o_string_stream.h\ + include/utils/perf_log.h src/solvers/diff_solver.$(obj-suffix):\ src/solvers/diff_solver.C\ include/base/auto_ptr.h\ @@ -8137,6 +8258,8 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ include/solvers/diff_system.h\ include/solvers/explicit_system.h\ include/solvers/implicit_system.h\ @@ -8192,6 +8315,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/solvers/eigen_system.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ @@ -8252,6 +8376,8 @@ include/numerics/sparse_matrix.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ include/solvers/diff_solver.h\ include/solvers/diff_system.h\ include/solvers/eigen_time_solver.h\ @@ -8305,6 +8431,7 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -8360,6 +8487,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -8417,6 +8545,8 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ include/solvers/diff_system.h\ include/solvers/euler2_solver.h\ include/solvers/explicit_system.h\ @@ -8471,6 +8601,8 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ include/solvers/diff_system.h\ include/solvers/euler_solver.h\ include/solvers/explicit_system.h\ @@ -8537,6 +8669,7 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/solvers/equation_systems.h\ include/solvers/exact_solution.h\ @@ -8589,12 +8722,90 @@ include/mesh/mesh_base.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/explicit_system.h\ include/solvers/system.h\ include/solvers/system_norm.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h +src/solvers/fem_context.$(obj-suffix):\ + src/solvers/fem_context.C\ + include/base/auto_ptr.h\ + include/base/dof_map.h\ + include/base/dof_object.h\ + include/base/libmesh.h\ + include/base/libmesh_base.h\ + include/base/libmesh_common.h\ + include/base/libmesh_config.h\ + include/base/libmesh_exceptions.h\ + include/base/libmesh_logging.h\ + include/base/multi_predicates.h\ + include/base/print_trace.h\ + include/base/reference_counted_object.h\ + include/base/reference_counter.h\ + include/base/single_predicates.h\ + include/base/variant_filter_iterator.h\ + include/enums/enum_elem_quality.h\ + include/enums/enum_elem_type.h\ + include/enums/enum_fe_family.h\ + include/enums/enum_inf_map_type.h\ + include/enums/enum_io_package.h\ + include/enums/enum_norm_type.h\ + include/enums/enum_order.h\ + include/enums/enum_quadrature_type.h\ + include/enums/enum_solver_package.h\ + include/enums/enum_xdr_mode.h\ + include/fe/fe_base.h\ + include/fe/fe_interface.h\ + include/fe/fe_type.h\ + include/geom/elem.h\ + include/geom/elem_range.h\ + include/geom/node.h\ + include/geom/point.h\ + include/geom/stored_range.h\ + include/mesh/mesh.h\ + include/mesh/mesh_base.h\ + include/mesh/parallel_mesh.h\ + include/mesh/serial_mesh.h\ + include/mesh/unstructured_mesh.h\ + include/numerics/dense_matrix.h\ + include/numerics/dense_matrix_base.h\ + include/numerics/dense_submatrix.h\ + include/numerics/dense_subvector.h\ + include/numerics/dense_vector.h\ + include/numerics/dense_vector_base.h\ + include/numerics/numeric_vector.h\ + include/numerics/sparse_matrix.h\ + include/numerics/tensor_value.h\ + include/numerics/type_tensor.h\ + include/numerics/type_vector.h\ + include/numerics/vector_value.h\ + include/parallel/threads.h\ + include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ + include/quadrature/quadrature.h\ + include/solvers/diff_context.h\ + include/solvers/diff_system.h\ + include/solvers/equation_systems.h\ + include/solvers/explicit_system.h\ + include/solvers/fem_context.h\ + include/solvers/fem_system.h\ + include/solvers/frequency_system.h\ + include/solvers/implicit_system.h\ + include/solvers/linear_implicit_system.h\ + include/solvers/newmark_system.h\ + include/solvers/steady_system.h\ + include/solvers/system.h\ + include/solvers/system_norm.h\ + include/solvers/time_solver.h\ + include/solvers/transient_system.h\ + include/solvers/unsteady_solver.h\ + include/utils/compare_types.h\ + include/utils/mapvector.h\ + include/utils/o_string_stream.h\ + include/utils/parameters.h\ + include/utils/perf_log.h src/solvers/fem_system.$(obj-suffix):\ src/solvers/fem_system.C\ include/base/auto_ptr.h\ @@ -8649,10 +8860,13 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ + include/solvers/diff_context.h\ include/solvers/diff_system.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ + include/solvers/fem_context.h\ include/solvers/fem_system.h\ include/solvers/frequency_system.h\ include/solvers/implicit_system.h\ @@ -8707,6 +8921,7 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -8764,6 +8979,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/solvers/explicit_system.h\ include/solvers/implicit_system.h\ include/solvers/system.h\ @@ -8809,6 +9025,7 @@ include/numerics/linear_solver.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -8859,6 +9076,7 @@ include/numerics/sparse_matrix.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -8921,6 +9139,8 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ include/solvers/diff_solver.h\ include/solvers/diff_system.h\ include/solvers/equation_systems.h\ @@ -8973,6 +9193,7 @@ include/numerics/nonlinear_solver.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -9037,6 +9258,8 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ include/solvers/diff_solver.h\ include/solvers/diff_system.h\ include/solvers/explicit_system.h\ @@ -9090,6 +9313,8 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ include/solvers/diff_solver.h\ include/solvers/diff_system.h\ include/solvers/explicit_system.h\ @@ -9148,6 +9373,7 @@ include/parallel/parallel.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ @@ -9206,6 +9432,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/system.h\ include/solvers/system_norm.h\ include/utils/compare_types.h\ @@ -9264,6 +9491,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/quadrature/quadrature.h\ include/quadrature/quadrature_gauss.h\ include/solvers/system.h\ @@ -9316,6 +9544,8 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ include/solvers/diff_solver.h\ include/solvers/diff_system.h\ include/solvers/explicit_system.h\ @@ -9366,6 +9596,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/solvers/explicit_system.h\ include/solvers/implicit_system.h\ include/solvers/linear_implicit_system.h\ @@ -9418,7 +9649,9 @@ include/numerics/numeric_vector.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/solvers/adaptive_time_solver.h\ + include/solvers/diff_context.h\ include/solvers/diff_system.h\ include/solvers/euler_solver.h\ include/solvers/explicit_system.h\ @@ -9477,6 +9710,8 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ + include/solvers/diff_context.h\ include/solvers/diff_solver.h\ include/solvers/diff_system.h\ include/solvers/explicit_system.h\ @@ -9537,6 +9772,7 @@ include/numerics/vector_value.h\ include/parallel/threads.h\ include/parallel/threads_allocators.h\ + include/partitioning/partitioner.h\ include/solvers/equation_systems.h\ include/solvers/explicit_system.h\ include/solvers/frequency_system.h\ @@ -9582,6 +9818,7 @@ include/numerics/type_vector.h\ include/parallel/parallel.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/location_maps.h\ include/utils/o_string_stream.h\ @@ -9680,6 +9917,7 @@ include/mesh/mesh_base.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h\ @@ -9714,6 +9952,7 @@ include/mesh/mesh_tools.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h\ @@ -9771,6 +10010,7 @@ include/mesh/mesh_tools.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h\ @@ -9804,6 +10044,7 @@ include/mesh/mesh_base.h\ include/numerics/type_vector.h\ include/parallel/threads.h\ + include/partitioning/partitioner.h\ include/utils/compare_types.h\ include/utils/o_string_stream.h\ include/utils/perf_log.h\ Index: examples/ex18/naviersystem.h =================================================================== --- examples/ex18/naviersystem.h (revision 3206) +++ examples/ex18/naviersystem.h (working copy) @@ -35,26 +35,27 @@ // System initialization virtual void init_data (); + // Context initialization + virtual void init_context(DiffContext &context); + // Element residual and jacobian calculations // Time dependent parts - virtual bool element_time_derivative (bool request_jacobian); + virtual bool element_time_derivative (bool request_jacobian, + DiffContext& context); // Constraint parts - virtual bool element_constraint (bool request_jacobian); - virtual bool side_constraint (bool request_jacobian); + virtual bool element_constraint (bool request_jacobian, + DiffContext& context); + virtual bool side_constraint (bool request_jacobian, + DiffContext& context); // Mass matrix part - virtual bool mass_residual (bool request_jacobian); + virtual bool mass_residual (bool request_jacobian, + DiffContext& context); // Indices for each variable; unsigned int p_var, u_var, v_var, w_var; - // Finite elements for the velocity and pressure on element interiors - FEBase *fe_velocity, *fe_pressure; - - // Finite element for the velocity on element sides - FEBase *fe_side_vel; - // The Reynolds number to solve for Real Reynolds; Index: examples/ex18/naviersystem.C =================================================================== --- examples/ex18/naviersystem.C (revision 3206) +++ examples/ex18/naviersystem.C (working copy) @@ -66,23 +66,6 @@ if (dim == 3) this->time_evolving(w_var); - // Get references to the finite elements we need - fe_velocity = element_fe[this->variable_type(u_var)]; - fe_pressure = element_fe[this->variable_type(p_var)]; - fe_side_vel = side_fe[this->variable_type(u_var)]; - - // To enable FE optimizations, we should prerequest all the data - // we will need to build the linear system. - fe_velocity->get_JxW(); - fe_velocity->get_phi(); - fe_velocity->get_dphi(); - fe_velocity->get_xyz(); - - fe_pressure->get_phi(); - - fe_side_vel->get_JxW(); - fe_side_vel->get_phi(); - // Check the input file for Reynolds number GetPot infile("navier.in"); Reynolds = infile("Reynolds", 1.); @@ -96,51 +79,77 @@ } -bool NavierSystem::element_time_derivative (bool request_jacobian) + +void NavierSystem::init_context(DiffContext &context) { + FEMContext &c = libmesh_cast_ref(context); + + // We should prerequest all the data + // we will need to build the linear system. + c.element_fe_var[u_var]->get_JxW(); + c.element_fe_var[u_var]->get_phi(); + c.element_fe_var[u_var]->get_dphi(); + c.element_fe_var[u_var]->get_xyz(); + + c.element_fe_var[p_var]->get_phi(); + + c.side_fe_var[u_var]->get_JxW(); + c.side_fe_var[u_var]->get_phi(); +} + + +bool NavierSystem::element_time_derivative (bool request_jacobian, + DiffContext &context) +{ + FEMContext &c = libmesh_cast_ref(context); + // First we get some references to cell-specific data that // will be used to assemble the linear system. // Element Jacobian * quadrature weights for interior integration - const std::vector &JxW = fe_velocity->get_JxW(); + const std::vector &JxW = + c.element_fe_var[u_var]->get_JxW(); // The velocity shape functions at interior quadrature points. - const std::vector >& phi = fe_velocity->get_phi(); + const std::vector >& phi = + c.element_fe_var[u_var]->get_phi(); // The velocity shape function gradients at interior // quadrature points. const std::vector >& dphi = - fe_velocity->get_dphi(); + c.element_fe_var[u_var]->get_dphi(); // The pressure shape functions at interior // quadrature points. - const std::vector >& psi = fe_pressure->get_phi(); + const std::vector >& psi = + c.element_fe_var[p_var]->get_phi(); // Physical location of the quadrature points - const std::vector& qpoint = fe_velocity->get_xyz(); + const std::vector& qpoint = + c.element_fe_var[u_var]->get_xyz(); // The number of local degrees of freedom in each variable - const unsigned int n_p_dofs = dof_indices_var[p_var].size(); - const unsigned int n_u_dofs = dof_indices_var[u_var].size(); - libmesh_assert (n_u_dofs == dof_indices_var[v_var].size()); + const unsigned int n_p_dofs = c.dof_indices_var[p_var].size(); + const unsigned int n_u_dofs = c.dof_indices_var[u_var].size(); + libmesh_assert (n_u_dofs == c.dof_indices_var[v_var].size()); // The subvectors and submatrices we need to fill: const unsigned int dim = this->get_mesh().mesh_dimension(); - DenseSubMatrix &Kuu = *elem_subjacobians[u_var][u_var]; - DenseSubMatrix &Kvv = *elem_subjacobians[v_var][v_var]; - DenseSubMatrix &Kww = *elem_subjacobians[w_var][w_var]; - DenseSubMatrix &Kuv = *elem_subjacobians[u_var][v_var]; - DenseSubMatrix &Kuw = *elem_subjacobians[u_var][w_var]; - DenseSubMatrix &Kvu = *elem_subjacobians[v_var][u_var]; - DenseSubMatrix &Kvw = *elem_subjacobians[v_var][w_var]; - DenseSubMatrix &Kwu = *elem_subjacobians[w_var][u_var]; - DenseSubMatrix &Kwv = *elem_subjacobians[w_var][v_var]; - DenseSubMatrix &Kup = *elem_subjacobians[u_var][p_var]; - DenseSubMatrix &Kvp = *elem_subjacobians[v_var][p_var]; - DenseSubMatrix &Kwp = *elem_subjacobians[w_var][p_var]; - DenseSubVector &Fu = *elem_subresiduals[u_var]; - DenseSubVector &Fv = *elem_subresiduals[v_var]; - DenseSubVector &Fw = *elem_subresiduals[w_var]; + DenseSubMatrix &Kuu = *c.elem_subjacobians[u_var][u_var]; + DenseSubMatrix &Kvv = *c.elem_subjacobians[v_var][v_var]; + DenseSubMatrix &Kww = *c.elem_subjacobians[w_var][w_var]; + DenseSubMatrix &Kuv = *c.elem_subjacobians[u_var][v_var]; + DenseSubMatrix &Kuw = *c.elem_subjacobians[u_var][w_var]; + DenseSubMatrix &Kvu = *c.elem_subjacobians[v_var][u_var]; + DenseSubMatrix &Kvw = *c.elem_subjacobians[v_var][w_var]; + DenseSubMatrix &Kwu = *c.elem_subjacobians[w_var][u_var]; + DenseSubMatrix &Kwv = *c.elem_subjacobians[w_var][v_var]; + DenseSubMatrix &Kup = *c.elem_subjacobians[u_var][p_var]; + DenseSubMatrix &Kvp = *c.elem_subjacobians[v_var][p_var]; + DenseSubMatrix &Kwp = *c.elem_subjacobians[w_var][p_var]; + DenseSubVector &Fu = *c.elem_subresiduals[u_var]; + DenseSubVector &Fv = *c.elem_subresiduals[v_var]; + DenseSubVector &Fw = *c.elem_subresiduals[w_var]; // Now we will build the element Jacobian and residual. // Constructing the residual requires the solution and its @@ -148,18 +157,18 @@ // calculated at each quadrature point by summing the // solution degree-of-freedom values by the appropriate // weight functions. - unsigned int n_qpoints = element_qrule->n_points(); + unsigned int n_qpoints = c.element_qrule->n_points(); for (unsigned int qp=0; qp != n_qpoints; qp++) { // Compute the solution & its gradient at the old Newton iterate - Number p = interior_value(p_var, qp), - u = interior_value(u_var, qp), - v = interior_value(v_var, qp), - w = interior_value(w_var, qp); - Gradient grad_u = interior_gradient(u_var, qp), - grad_v = interior_gradient(v_var, qp), - grad_w = interior_gradient(w_var, qp); + Number p = c.interior_value(p_var, qp), + u = c.interior_value(u_var, qp), + v = c.interior_value(v_var, qp), + w = c.interior_value(w_var, qp); + Gradient grad_u = c.interior_gradient(u_var, qp), + grad_v = c.interior_gradient(v_var, qp), + grad_w = c.interior_gradient(w_var, qp); // Definitions for convenience. It is sometimes simpler to do a // dot product if you have the full vector at your disposal. @@ -265,42 +274,46 @@ -bool NavierSystem::element_constraint (bool request_jacobian) +bool NavierSystem::element_constraint (bool request_jacobian, + DiffContext &context) { + FEMContext &c = libmesh_cast_ref(context); + // Here we define some references to cell-specific data that // will be used to assemble the linear system. // Element Jacobian * quadrature weight for interior integration - const std::vector &JxW = fe_velocity->get_JxW(); + const std::vector &JxW = c.element_fe_var[u_var]->get_JxW(); // The velocity shape function gradients at interior // quadrature points. const std::vector >& dphi = - fe_velocity->get_dphi(); + c.element_fe_var[u_var]->get_dphi(); // The pressure shape functions at interior // quadrature points. - const std::vector >& psi = fe_pressure->get_phi(); + const std::vector >& psi = + c.element_fe_var[p_var]->get_phi(); // The number of local degrees of freedom in each variable - const unsigned int n_u_dofs = dof_indices_var[u_var].size(); - const unsigned int n_p_dofs = dof_indices_var[p_var].size(); + const unsigned int n_u_dofs = c.dof_indices_var[u_var].size(); + const unsigned int n_p_dofs = c.dof_indices_var[p_var].size(); // The subvectors and submatrices we need to fill: const unsigned int dim = this->get_mesh().mesh_dimension(); - DenseSubMatrix &Kpu = *elem_subjacobians[p_var][u_var]; - DenseSubMatrix &Kpv = *elem_subjacobians[p_var][v_var]; - DenseSubMatrix &Kpw = *elem_subjacobians[p_var][w_var]; - DenseSubVector &Fp = *elem_subresiduals[p_var]; + DenseSubMatrix &Kpu = *c.elem_subjacobians[p_var][u_var]; + DenseSubMatrix &Kpv = *c.elem_subjacobians[p_var][v_var]; + DenseSubMatrix &Kpw = *c.elem_subjacobians[p_var][w_var]; + DenseSubVector &Fp = *c.elem_subresiduals[p_var]; // Add the constraint given by the continuity equation - unsigned int n_qpoints = element_qrule->n_points(); + unsigned int n_qpoints = c.element_qrule->n_points(); for (unsigned int qp=0; qp != n_qpoints; qp++) { // Compute the velocity gradient at the old Newton iterate - Gradient grad_u = interior_gradient(u_var, qp), - grad_v = interior_gradient(v_var, qp), - grad_w = interior_gradient(w_var, qp); + Gradient grad_u = c.interior_gradient(u_var, qp), + grad_v = c.interior_gradient(v_var, qp), + grad_w = c.interior_gradient(w_var, qp); // Now a loop over the pressure degrees of freedom. This // computes the contributions of the continuity equation. @@ -328,30 +341,33 @@ -bool NavierSystem::side_constraint (bool request_jacobian) +bool NavierSystem::side_constraint (bool request_jacobian, + DiffContext &context) { + FEMContext &c = libmesh_cast_ref(context); + // Here we define some references to cell-specific data that // will be used to assemble the linear system. // Element Jacobian * quadrature weight for side integration - const std::vector &JxW_side = fe_side_vel->get_JxW(); + const std::vector &JxW_side = c.side_fe_var[u_var]->get_JxW(); // The velocity shape functions at side quadrature points. const std::vector >& phi_side = - fe_side_vel->get_phi(); + c.side_fe_var[u_var]->get_phi(); // The number of local degrees of freedom in u and v - const unsigned int n_u_dofs = dof_indices_var[u_var].size(); + const unsigned int n_u_dofs = c.dof_indices_var[u_var].size(); // The subvectors and submatrices we need to fill: const unsigned int dim = this->get_mesh().mesh_dimension(); - DenseSubMatrix &Kuu = *elem_subjacobians[u_var][u_var]; - DenseSubMatrix &Kvv = *elem_subjacobians[v_var][v_var]; - DenseSubMatrix &Kww = *elem_subjacobians[w_var][w_var]; + DenseSubMatrix &Kuu = *c.elem_subjacobians[u_var][u_var]; + DenseSubMatrix &Kvv = *c.elem_subjacobians[v_var][v_var]; + DenseSubMatrix &Kww = *c.elem_subjacobians[w_var][w_var]; - DenseSubVector &Fu = *elem_subresiduals[u_var]; - DenseSubVector &Fv = *elem_subresiduals[v_var]; - DenseSubVector &Fw = *elem_subresiduals[w_var]; + DenseSubVector &Fu = *c.elem_subresiduals[u_var]; + DenseSubVector &Fv = *c.elem_subresiduals[v_var]; + DenseSubVector &Fw = *c.elem_subresiduals[w_var]; // For this example we will use Dirichlet velocity boundary // conditions imposed at each timestep via the penalty method. @@ -359,19 +375,19 @@ // The penalty value. \f$ \frac{1}{\epsilon} \f$ const Real penalty = 1.e10; - unsigned int n_sidepoints = side_qrule->n_points(); + unsigned int n_sidepoints = c.side_qrule->n_points(); for (unsigned int qp=0; qp != n_sidepoints; qp++) { // Compute the solution at the old Newton iterate - Number u = side_value(u_var, qp), - v = side_value(v_var, qp), - w = side_value(w_var, qp); + Number u = c.side_value(u_var, qp), + v = c.side_value(v_var, qp), + w = c.side_value(w_var, qp); // Set u = 1 on the top boundary, (which build_square() has // given boundary id 2, and build_cube() has called 5), // u = 0 everywhere else short int boundary_id = - this->get_mesh().boundary_info->boundary_id(elem, side); + this->get_mesh().boundary_info->boundary_id(c.elem, c.side); libmesh_assert (boundary_id != BoundaryInfo::invalid_id); const short int top_id = (dim==3) ? 5 : 2; @@ -410,27 +426,27 @@ } // Pin p = 0 at the origin - if (elem->contains_point(Point(0.,0.))) + if (c.elem->contains_point(Point(0.,0.))) { // The pressure penalty value. \f$ \frac{1}{\epsilon} \f$ const Real penalty = 1.e9; - DenseSubMatrix &Kpp = *elem_subjacobians[p_var][p_var]; - DenseSubVector &Fp = *elem_subresiduals[p_var]; - const unsigned int n_p_dofs = dof_indices_var[p_var].size(); + DenseSubMatrix &Kpp = *c.elem_subjacobians[p_var][p_var]; + DenseSubVector &Fp = *c.elem_subresiduals[p_var]; + const unsigned int n_p_dofs = c.dof_indices_var[p_var].size(); Point zero(0.); - Number p = point_value(p_var, zero); + Number p = c.point_value(p_var, zero); Number p_value = 0.; unsigned int dim = get_mesh().mesh_dimension(); - FEType fe_type = element_fe_var[p_var]->get_fe_type(); - Point p_master = FEInterface::inverse_map(dim, fe_type, elem, zero); + FEType fe_type = c.element_fe_var[p_var]->get_fe_type(); + Point p_master = FEInterface::inverse_map(dim, fe_type, c.elem, zero); std::vector point_phi(n_p_dofs); for (unsigned int i=0; i != n_p_dofs; i++) { - point_phi[i] = FEInterface::shape(dim, fe_type, elem, i, p_master); + point_phi[i] = FEInterface::shape(dim, fe_type, c.elem, i, p_master); } for (unsigned int i=0; i != n_p_dofs; i++) @@ -455,35 +471,39 @@ // Reynolds number (or choose a more complicated non-dimensionalization // of time), but this gives us an opportunity to demonstrate overriding // FEMSystem::mass_residual() -bool NavierSystem::mass_residual (bool request_jacobian) +bool NavierSystem::mass_residual (bool request_jacobian, + DiffContext &context) { + FEMContext &c = libmesh_cast_ref(context); + // The subvectors and submatrices we need to fill: const unsigned int dim = this->get_mesh().mesh_dimension(); // Element Jacobian * quadrature weight for interior integration - const std::vector &JxW = fe_velocity->get_JxW(); + const std::vector &JxW = c.element_fe_var[u_var]->get_JxW(); // The velocity shape functions at interior quadrature points. - const std::vector >& phi = fe_velocity->get_phi(); + const std::vector >& phi = + c.element_fe_var[u_var]->get_phi(); // The subvectors and submatrices we need to fill: - DenseSubVector &Fu = *elem_subresiduals[u_var]; - DenseSubVector &Fv = *elem_subresiduals[v_var]; - DenseSubVector &Fw = *elem_subresiduals[w_var]; - DenseSubMatrix &Kuu = *elem_subjacobians[u_var][u_var]; - DenseSubMatrix &Kvv = *elem_subjacobians[v_var][v_var]; - DenseSubMatrix &Kww = *elem_subjacobians[w_var][w_var]; + DenseSubVector &Fu = *c.elem_subresiduals[u_var]; + DenseSubVector &Fv = *c.elem_subresiduals[v_var]; + DenseSubVector &Fw = *c.elem_subresiduals[w_var]; + DenseSubMatrix &Kuu = *c.elem_subjacobians[u_var][u_var]; + DenseSubMatrix &Kvv = *c.elem_subjacobians[v_var][v_var]; + DenseSubMatrix &Kww = *c.elem_subjacobians[w_var][w_var]; // The number of local degrees of freedom in velocity - const unsigned int n_u_dofs = dof_indices_var[u_var].size(); + const unsigned int n_u_dofs = c.dof_indices_var[u_var].size(); - unsigned int n_qpoints = element_qrule->n_points(); + unsigned int n_qpoints = c.element_qrule->n_points(); for (unsigned int qp = 0; qp != n_qpoints; ++qp) { - Number u = interior_value(u_var, qp), - v = interior_value(v_var, qp), - w = interior_value(w_var, qp); + Number u = c.interior_value(u_var, qp), + v = c.interior_value(v_var, qp), + w = c.interior_value(w_var, qp); // We pull as many calculations as possible outside of loops Number JxWxRe = JxW[qp] * Reynolds; @@ -498,9 +518,9 @@ if (dim == 3) Fw(i) += JxWxRexW * phi[i][qp]; - if (request_jacobian && elem_solution_derivative) + if (request_jacobian && c.elem_solution_derivative) { - libmesh_assert (elem_solution_derivative == 1.0); + libmesh_assert (c.elem_solution_derivative == 1.0); Number JxWxRexPhiI = JxWxRe * phi[i][qp]; Number JxWxRexPhiII = JxWxRexPhiI * phi[i][qp];