Index: include/systems/fem_system.h =================================================================== --- include/systems/fem_system.h (revision 6309) +++ include/systems/fem_system.h (working copy) @@ -23,6 +23,7 @@ // Local Includes #include "libmesh/diff_system.h" +#include "libmesh/fem_physics.h" // C++ includes #include @@ -52,7 +53,8 @@ // ------------------------------------------------------------ // FEMSystem class definition -class FEMSystem : public DifferentiableSystem +class FEMSystem : public DifferentiableSystem, + public FEMPhysics { public: @@ -103,72 +105,6 @@ virtual void solve (); /** - * Tells the FEMSystem that variable \p var is evolving with - * respect to time. In general, the user's init() function - * should call time_evolving() for any variables which - * behave like du/dt = F(u), and should not call time_evolving() - * for any variables which behave like 0 = G(u). - * - * Most derived systems will not have to reimplment this function; however - * any system which reimplements mass_residual() may have to reimplement - * time_evolving() to prepare data structures. - */ - virtual void time_evolving (unsigned int var); - - /** - * Tells the FEMSystem that system \p sys contains the - * isoparametric Lagrangian variables which correspond to the - * coordinates of mesh nodes, in problems where the mesh itself is - * expected to move in time. - * - * The system with mesh coordinate data (which may be \p this system - * itself, for fully coupled moving mesh problems) is currently - * assumed to have new (end of time step) mesh coordinates stored in - * solution, old (beginning of time step) mesh coordinates stored in - * _old_nonlinear_solution, and constant velocity motion during each - * time step. - * - * Activating this function ensures that local (but not neighbor!) element - * geometry is correctly repositioned when evaluating element residuals. - * - * Currently \p sys must be \p *this for a tightly coupled moving - * mesh problem or NULL to stop mesh movement; loosely coupled - * moving mesh problems are not implemented. - * - * This code is experimental. "Trust but verify, and not in that - * order" - */ - virtual void set_mesh_system(System* sys); - - /** - * Tells the FEMSystem that variable \p var from the mesh system - * should be used to update the x coordinate of mesh nodes, in problems where - * the mesh itself is expected to move in time. - * - * The system with mesh coordinate data (which may be this system itself, for - * fully coupled moving mesh problems) is currently assumed to have new (end - * of time step) mesh coordinates stored in solution, old (beginning of time - * step) mesh coordinates stored in _old_nonlinear_solution, and constant - * velocity motion during each time step. - * - * Activating this function ensures that local (but not neighbor!) element - * geometry is correctly repositioned when evaluating element residuals. - */ - virtual void set_mesh_x_var(unsigned int var); - - /** - * Tells the FEMSystem that variable \p var from the mesh system - * should be used to update the y coordinate of mesh nodes. - */ - virtual void set_mesh_y_var(unsigned int var); - - /** - * Tells the FEMSystem that variable \p var from the mesh system - * should be used to update the z coordinate of mesh nodes. - */ - virtual void set_mesh_z_var(unsigned int var); - - /** * Tells the FEMSystem to set the degree of freedom coefficients * which should correspond to mesh nodal coordinates. */ @@ -181,39 +117,6 @@ void mesh_position_set(); /** - * Adds a pseudo-convection contribution on \p elem to - * elem_residual, if the nodes of \p elem are being translated by a - * moving mesh. - * - * This function assumes that the user's time derivative equations - * (except for any equations involving unknown mesh xyz coordinates - * themselves) are expressed in an Eulerian frame of reference, and - * that the user is satisfied with an unstabilized convection term. - * Lagrangian equations will probably require overriding - * eulerian_residual() with a blank function; ALE or stabilized - * formulations will require reimplementing eulerian_residual() - * entirely. - */ - virtual bool eulerian_residual (bool request_jacobian, - DiffContext &context); - - /** - * Adds a mass vector contribution on \p elem to elem_residual. - * If this method receives request_jacobian = true, then it - * should compute elem_jacobian and return true if possible. If - * elem_jacobian has not been computed then the method should - * return false. - * - * Most problems can use the FEMSystem::mass_residual implementation, - * which calculates the residual (u, phi_i) and jacobian (phi_i, phi_j); - * few users will need to reimplement this themselves. Using a custom - * mass matrix (e.g. for divergence-free elements or mass lumping) - * requires reimplementing mass_residual(). - */ - virtual bool mass_residual (bool request_jacobian, - DiffContext &context); - - /** * Builds a FEMContext object with enough information to do * evaluations on each element. * @@ -323,54 +226,6 @@ }; - -// ------------------------------------------------------------ -// FEMSystem inline methods - - - -inline -void FEMSystem::set_mesh_system(System* sys) -{ - // For now we assume that we're doing fully coupled mesh motion - if (sys && sys != this) - libmesh_not_implemented(); - - // For the foreseeable future we'll assume that we keep these - // Systems in the same EquationSystems - libmesh_assert_equal_to (&this->get_equation_systems(), - &sys->get_equation_systems()); - - // And for the immediate future this code may not even work - libmesh_experimental(); - - _mesh_sys = sys; -} - - - -inline -void FEMSystem::set_mesh_x_var (unsigned int var) -{ - _mesh_x_var = var; -} - - - -inline -void FEMSystem::set_mesh_y_var (unsigned int var) -{ - _mesh_y_var = var; -} - - - -inline -void FEMSystem::set_mesh_z_var (unsigned int var) -{ - _mesh_z_var = var; -} - } // namespace libMesh Index: include/systems/diff_system.h =================================================================== --- include/systems/diff_system.h (revision 6309) +++ include/systems/diff_system.h (working copy) @@ -54,8 +54,8 @@ // DifferentiableSystem class definition class DifferentiableSystem : public ImplicitSystem, - public DifferentiablePhysics, - public DifferentiableQoI + public virtual DifferentiablePhysics, + public virtual DifferentiableQoI { public: @@ -134,14 +134,37 @@ virtual void solve (); /** - * Force the user to override clone for DifferentiableQoI + * We don't allow systems to be attached to each other */ + virtual AutoPtr clone_physics() + { libmesh_error(); + // dummy + return AutoPtr(this); } + + /** + * We don't allow systems to be attached to each other + */ virtual AutoPtr clone() { libmesh_error(); // dummy return AutoPtr(this); } /** + * Returns const reference to DifferentiablePhysics object. Note + * that if no external Physics object is attached, the default is + * this. + */ + const DifferentiablePhysics* get_physics() + { return this->_diff_physics; } + + /** + * Attach external Physics object. + */ + void attach_physics( DifferentiablePhysics* physics_in ) + { this->_diff_physics = (physics_in->clone_physics()).release(); + this->_diff_physics->init_physics(*this);} + + /** * Returns const reference to DifferentiableQoI object. Note that if no external * QoI object is attached, the default is this. */ @@ -238,15 +261,15 @@ */ bool print_element_jacobians; +protected: + /** * Pointer to object to use for physics assembly evaluations. * Defaults to \p this for backwards compatibility; in the future * users should create separate physics objects. */ - DifferentiablePhysics *diff_physics; + DifferentiablePhysics *_diff_physics; -protected: - /** * Pointer to object to use for quantity of interest assembly * evaluations. Defaults to \p this for backwards compatibility; in Index: include/systems/diff_context.h =================================================================== --- include/systems/diff_context.h (revision 6309) +++ include/systems/diff_context.h (working copy) @@ -86,6 +86,11 @@ virtual void elem_edge_reinit(Real) {} /** + * Number of variables in solution. + */ + unsigned int n_vars() const { return dof_indices_var.size(); } + + /** * Accessor for element solution. */ const DenseVector& get_elem_solution() const Index: include/systems/system.h =================================================================== --- include/systems/system.h (revision 6309) +++ include/systems/system.h (working copy) @@ -1238,31 +1238,6 @@ virtual void prolong_vectors (); /** - * Returns a reference to the system with variables corresponding to - * mesh nodal coordinates, or NULL if the mesh is fixed. - * Useful for ALE calculations. - */ - const System* get_mesh_system() const; - - /** - * Returns the variable number corresponding to the - * mesh x coordinate. Useful for ALE calculations. - */ - unsigned int get_mesh_x_var() const; - - /** - * Returns the variable number corresponding to the - * mesh y coordinate. Useful for ALE calculations. - */ - unsigned int get_mesh_y_var() const; - - /** - * Returns the variable number corresponding to the - * mesh z coordinate. Useful for ALE calculations. - */ - unsigned int get_mesh_z_var() const; - - /** * Flag which tells the system to whether or not to * call the user assembly function during each call to solve(). * By default, every call to solve() begins with a call to the @@ -1456,16 +1431,6 @@ void project_vector (const NumericVector&, NumericVector&) const; - /** - * System from which to acquire moving mesh information - */ - System *_mesh_sys; - - /** - * Variables from which to acquire moving mesh information - */ - unsigned int _mesh_x_var, _mesh_y_var, _mesh_z_var; - private: /** * This isn't a copyable object, so let's make sure nobody tries. @@ -1995,30 +1960,7 @@ libmesh_not_implemented(); } -inline -const System* System::get_mesh_system() const -{ - return _mesh_sys; -} -inline -unsigned int System::get_mesh_x_var() const -{ - return _mesh_x_var; -} - -inline -unsigned int System::get_mesh_y_var() const -{ - return _mesh_y_var; -} - -inline -unsigned int System::get_mesh_z_var() const -{ - return _mesh_z_var; -} - } // namespace libMesh #endif // #define __system_h__ Index: include/physics/fem_physics.h =================================================================== --- include/physics/fem_physics.h (revision 6309) +++ include/physics/fem_physics.h (working copy) @@ -18,25 +18,17 @@ -#ifndef __diff_physics_h__ -#define __diff_physics_h__ +#ifndef __fem_physics_h__ +#define __fem_physics_h__ // Local Includes -#include "libmesh/auto_ptr.h" -#include "libmesh/diff_context.h" +#include "libmesh/diff_physics.h" // C++ includes -#include namespace libMesh { -// Forward Declarations -class System; -class TimeSolver; - -template class NumericVector; - /** * This class provides a specific system class. It aims * to generalize any system, linear or nonlinear, which @@ -46,141 +38,128 @@ * which is still experimental. Users of this framework should * beware of bugs and future API changes. * - * @author Roy H. Stogner 2006 + * @author Roy H. Stogner 2012 */ // ------------------------------------------------------------ -// DifferentiableSystem class definition +// Finite Element Method Physics class definition -class DifferentiablePhysics +class FEMPhysics : public virtual DifferentiablePhysics { public: /** - * Constructor. Optionally initializes required - * data structures. + * Constructor. */ - DifferentiablePhysics () : compute_internal_sides(false) {} + FEMPhysics () : + DifferentiablePhysics(), + _mesh_sys (NULL), + _mesh_x_var (libMesh::invalid_uint), + _mesh_y_var (libMesh::invalid_uint), + _mesh_z_var (libMesh::invalid_uint) + {} /** * Destructor. */ - virtual ~DifferentiablePhysics (); + virtual ~FEMPhysics () {} /** - * Clear any data structures associated with the physics. + * Tells the FEMPhysics that system \p sys contains the + * isoparametric Lagrangian variables which correspond to the + * coordinates of mesh nodes, in problems where the mesh itself is + * expected to move in time. + * + * The system with mesh coordinate data (which may be \p this system + * itself, for fully coupled moving mesh problems) is currently + * assumed to have new (end of time step) mesh coordinates stored in + * solution, old (beginning of time step) mesh coordinates stored in + * _old_nonlinear_solution, and constant velocity motion during each + * time step. + * + * Activating this function ensures that local (but not neighbor!) element + * geometry is correctly repositioned when evaluating element residuals. + * + * Currently \p sys must be \p *this for a tightly coupled moving + * mesh problem or NULL to stop mesh movement; loosely coupled + * moving mesh problems are not implemented. + * + * This code is experimental. "Trust but verify, and not in that + * order" */ - virtual void clear_physics (); + virtual void set_mesh_system(System* sys); /** - * Initialize any data structures associated with the physics. + * Returns a reference to the system with variables corresponding to + * mesh nodal coordinates, or NULL if the mesh is fixed. + * Useful for ALE calculations. */ - virtual void init_physics (const System& sys); + const System* get_mesh_system() const; /** - * Adds the time derivative contribution on \p elem to elem_residual. - * If this method receives request_jacobian = true, then it - * should compute elem_jacobian and return true if possible. If - * elem_jacobian has not been computed then the method should - * return false. + * Tells the FEMPhysics that variable \p var from the mesh system + * should be used to update the x coordinate of mesh nodes, in problems where + * the mesh itself is expected to move in time. * - * Users need to reimplement this for their particular PDE. - * To implement the physics model du/dt = F(u), the user should - * examine u = elem_solution and add (F(u), phi_i) to - * elem_residual. + * The system with mesh coordinate data (which may be this system itself, for + * fully coupled moving mesh problems) is currently assumed to have new (end + * of time step) mesh coordinates stored in solution, old (beginning of time + * step) mesh coordinates stored in _old_nonlinear_solution, and constant + * velocity motion during each time step. + * + * Activating this function ensures that local (but not neighbor!) element + * geometry is correctly repositioned when evaluating element residuals. */ - virtual bool element_time_derivative (bool request_jacobian, - DiffContext &) { - return request_jacobian; - } + virtual void set_mesh_x_var(unsigned int var); /** - * Adds the constraint contribution on \p elem to elem_residual. - * If this method receives request_jacobian = true, then it - * should compute elem_jacobian and return true if possible. If - * elem_jacobian has not been computed then the method should - * return false. - * - * Users may need to reimplement this for their particular PDE. - * To implement the constraint 0 = G(u), the user should - * examine u = elem_solution and add (G(u), phi_i) to - * elem_residual. + * Returns the variable number corresponding to the + * mesh x coordinate. Useful for ALE calculations. */ - virtual bool element_constraint (bool request_jacobian, - DiffContext &) { - return request_jacobian; - } + unsigned int get_mesh_x_var() const; /** - * \p compute_internal_sides is false by default, indicating that - * side_* computations will only be done on boundary sides. If - * compute_internal_sides is true, computations will be done - * on sides between elements as well. + * Tells the FEMPhysics that variable \p var from the mesh system + * should be used to update the y coordinate of mesh nodes. */ - bool compute_internal_sides; + virtual void set_mesh_y_var(unsigned int var); /** - * Adds the time derivative contribution on \p side of \p elem to - * elem_residual. - * If this method receives request_jacobian = true, then it - * should compute elem_jacobian and return true if possible. If - * elem_jacobian has not been computed then the method should - * return false. - * - * Users may need to reimplement this for their particular PDE - * depending on the boundary conditions. + * Returns the variable number corresponding to the + * mesh y coordinate. Useful for ALE calculations. */ - virtual bool side_time_derivative (bool request_jacobian, - DiffContext &) { - return request_jacobian; - } + unsigned int get_mesh_y_var() const; /** - * Adds the time derivative contribution on \p side of \p elem to - * elem_residual. - * If this method receives request_jacobian = true, then it - * should compute elem_jacobian and return true if possible. If - * elem_jacobian has not been computed then the method should - * return false. - * - * Users may need to reimplement this for their particular PDE - * depending on the boundary conditions. + * Tells the FEMPhysics that variable \p var from the mesh system + * should be used to update the z coordinate of mesh nodes. */ - virtual bool side_constraint (bool request_jacobian, - DiffContext &) { - return request_jacobian; - } + virtual void set_mesh_z_var(unsigned int var); /** - * Tells the DiffSystem that variable var is evolving with - * respect to time. In general, the user's init() function - * should call time_evolving() for any variables which - * behave like du/dt = F(u), and should not call time_evolving() - * for any variables which behave like 0 = G(u). - * - * Most derived systems will not have to reimplment this function; however - * any system which reimplements mass_residual() may have to reimplement - * time_evolving() to prepare data structures. + * Returns the variable number corresponding to the + * mesh z coordinate. Useful for ALE calculations. */ - virtual void time_evolving (unsigned int var) { - if (_time_evolving.size() <= var) - _time_evolving.resize(var+1, false); - _time_evolving[var] = true; - } + unsigned int get_mesh_z_var() const; /** * Adds a pseudo-convection contribution on \p elem to * elem_residual, if the nodes of \p elem are being translated by a * moving mesh. * - * The library provides a basic implementation in - * FEMSystem::eulerian_residual() + * This function assumes that the user's time derivative equations + * (except for any equations involving unknown mesh xyz coordinates + * themselves) are expressed in an Eulerian frame of reference, and + * that the user is satisfied with an unstabilized convection term. + * Lagrangian equations will probably require overriding + * eulerian_residual() with a blank function; ALE or stabilized + * formulations will require reimplementing eulerian_residual() + * entirely. */ virtual bool eulerian_residual (bool request_jacobian, - DiffContext &) { - return request_jacobian; - } + DiffContext &context); + /** * Adds a mass vector contribution on \p elem to elem_residual. * If this method receives request_jacobian = true, then it @@ -189,48 +168,101 @@ * return false. * * Most problems can use the reimplementation in - * FEMSystem::mass_residual; few users will need to reimplement + * FEMPhysics::mass_residual; few users will need to reimplement * this themselves. */ virtual bool mass_residual (bool request_jacobian, - DiffContext &) { - return request_jacobian; - } + DiffContext &); +protected: + /** - * Adds a mass vector contribution on \p side of \p elem to - * elem_residual. - * If this method receives request_jacobian = true, then it - * should compute elem_jacobian and return true if possible. If - * elem_jacobian has not been computed then the method should - * return false. - * - * For most problems, the default implementation of "do nothing" - * is correct; users with boundary conditions including time - * derivatives may need to reimplement this themselves. + * System from which to acquire moving mesh information */ - virtual bool side_mass_residual (bool request_jacobian, - DiffContext &) { - return request_jacobian; - } + System *_mesh_sys; - /* - * 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 &) {} - -protected: - /** - * Stores bools to tell us which variables are evolving - * in time and which are just constraints + * Variables from which to acquire moving mesh information */ - std::vector _time_evolving; + unsigned int _mesh_x_var, _mesh_y_var, _mesh_z_var; }; + + +// ------------------------------------------------------------ +// FEMPhysics inline methods + + +inline +void FEMPhysics::set_mesh_system(System* sys) +{ + // For now we assume that we're doing fully coupled mesh motion +// if (sys && sys != this) +// libmesh_not_implemented(); + + // For the foreseeable future we'll assume that we keep these + // Systems in the same EquationSystems + // libmesh_assert_equal_to (&this->get_equation_systems(), + // &sys->get_equation_systems()); + + // And for the immediate future this code may not even work + libmesh_experimental(); + + _mesh_sys = sys; +} + + + +inline +void FEMPhysics::set_mesh_x_var (unsigned int var) +{ + _mesh_x_var = var; +} + + + +inline +void FEMPhysics::set_mesh_y_var (unsigned int var) +{ + _mesh_y_var = var; +} + + + +inline +void FEMPhysics::set_mesh_z_var (unsigned int var) +{ + _mesh_z_var = var; +} + + + +inline +const System* FEMPhysics::get_mesh_system() const +{ + return _mesh_sys; +} + +inline +unsigned int FEMPhysics::get_mesh_x_var() const +{ + return _mesh_x_var; +} + +inline +unsigned int FEMPhysics::get_mesh_y_var() const +{ + return _mesh_y_var; +} + +inline +unsigned int FEMPhysics::get_mesh_z_var() const +{ + return _mesh_z_var; +} + + + } // namespace libMesh Index: include/physics/diff_physics.h =================================================================== --- include/physics/diff_physics.h (revision 6309) +++ include/physics/diff_physics.h (working copy) @@ -68,6 +68,11 @@ virtual ~DifferentiablePhysics (); /** + * Copy of this object. User should override to copy any needed state. + */ + virtual AutoPtr clone_physics() = 0; + + /** * Clear any data structures associated with the physics. */ virtual void clear_physics (); @@ -169,6 +174,17 @@ } /** + * Returns true iff variable \p var is evolving with + * respect to time. In general, the user's init() function + * should have set time_evolving() for any variables which + * behave like du/dt = F(u), and should not call time_evolving() + * for any variables which behave like 0 = G(u). + */ + bool is_time_evolving (unsigned int var) const { + return _time_evolving[var]; + } + + /** * Adds a pseudo-convection contribution on \p elem to * elem_residual, if the nodes of \p elem are being translated by a * moving mesh. Index: include/libmesh/fem_physics.h =================================================================== --- include/libmesh/fem_physics.h (revision 0) +++ include/libmesh/fem_physics.h (revision 0) @@ -0,0 +1 @@ +link ../physics/fem_physics.h \ No newline at end of file Property changes on: include/libmesh/fem_physics.h ___________________________________________________________________ Added: svn:special + * Index: src/physics/fem_physics.C =================================================================== --- src/physics/fem_physics.C (revision 6309) +++ src/physics/fem_physics.C (working copy) @@ -33,967 +33,19 @@ #include "libmesh/unsteady_solver.h" // For eulerian_residual -namespace { - using namespace libMesh; - - // give this guy some scope since there - // is underlying vector allocation upon - // creation/deletion - ConstElemRange elem_range; - - typedef Threads::spin_mutex femsystem_mutex; - femsystem_mutex assembly_mutex; - - class AssemblyContributions - { - public: - /** - * constructor to set context - */ - AssemblyContributions(FEMSystem &sys, - bool get_residual, - bool get_jacobian) : - _sys(sys), - _get_residual(get_residual), - _get_jacobian(get_jacobian) {} - - /** - * operator() for use with Threads::parallel_for(). - */ - void operator()(const ConstElemRange &range) const - { - AutoPtr con = _sys.build_context(); - FEMContext &_femcontext = libmesh_cast_ref(*con); - _sys.init_context(_femcontext); - - for (ConstElemRange::const_iterator elem_it = range.begin(); - elem_it != range.end(); ++elem_it) - { - Elem *el = const_cast(*elem_it); - - _femcontext.pre_fe_reinit(_sys, el); - _femcontext.elem_fe_reinit(); - - bool jacobian_computed = - _sys.time_solver->element_residual(_get_jacobian, _femcontext); - - // 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_equal_to (_femcontext.elem_jacobian.l1_norm(), 0.0); - // Logging of numerical jacobians is done separately - _sys.numerical_elem_jacobian(_femcontext); - } - - // Compute a numeric jacobian if we're asked to verify the - // analytic jacobian we got - if (_get_jacobian && jacobian_computed && - _sys.verify_analytic_jacobians != 0.0) - { - DenseMatrix analytic_jacobian(_femcontext.elem_jacobian); - - _femcontext.elem_jacobian.zero(); - // Logging of numerical jacobians is done separately - _sys.numerical_elem_jacobian(_femcontext); - - Real analytic_norm = analytic_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(_femcontext.elem_jacobian); - - // The matrix "analytic_jacobian" will now hold the error matrix - analytic_jacobian.add(-1.0, _femcontext.elem_jacobian); - Real error_norm = analytic_jacobian.l1_norm(); - - Real relative_error = error_norm / - std::max(analytic_norm, numerical_norm); - - if (relative_error > _sys.verify_analytic_jacobians) - { - libMesh::err << "Relative error " << relative_error - << " detected in analytic jacobian on element " - << _femcontext.elem->id() << '!' << std::endl; - - unsigned int old_precision = libMesh::out.precision(); - libMesh::out.precision(16); - libMesh::out << "J_analytic " << _femcontext.elem->id() << " = " - << _femcontext.elem_jacobian << std::endl; - analytic_jacobian.add(1.0, _femcontext.elem_jacobian); - libMesh::out << "J_numeric " << _femcontext.elem->id() << " = " - << analytic_jacobian << std::endl; - - libMesh::out.precision(old_precision); - - libmesh_error(); - } - } - - for (_femcontext.side = 0; - _femcontext.side != _femcontext.elem->n_sides(); - ++_femcontext.side) - { - // Don't compute on non-boundary sides unless requested - if (!_sys.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 - _femcontext.side_fe_reinit(); - - DenseMatrix old_jacobian; - // If we're in DEBUG mode, we should always verify that the - // user's side_residual function doesn't alter our existing - // jacobian and then lie about it -#ifndef DEBUG - // Even if we're not in DEBUG mode, when we're verifying - // analytic jacobians we'll want to verify each side's - // jacobian contribution separately. - /* PB: We also need to account for the case when the user wants to - use numerical Jacobians and not analytic Jacobians */ - if ( (_sys.verify_analytic_jacobians != 0.0 && _get_jacobian) || - (!jacobian_computed && _get_jacobian) ) -#endif // ifndef DEBUG - { - old_jacobian = _femcontext.elem_jacobian; - _femcontext.elem_jacobian.zero(); - } - jacobian_computed = - _sys.time_solver->side_residual(_get_jacobian, _femcontext); - - // Compute a numeric jacobian if we have to - if (_get_jacobian && !jacobian_computed) - { - // In DEBUG mode, we've already set elem_jacobian == 0, - // so we can make sure side_residual didn't compute a - // jacobian and lie about it -#ifdef DEBUG - libmesh_assert_equal_to (_femcontext.elem_jacobian.l1_norm(), 0.0); -#endif - // Logging of numerical jacobians is done separately - _sys.numerical_side_jacobian(_femcontext); - - // Add back in element interior numerical Jacobian - _femcontext.elem_jacobian += old_jacobian; - } - - // Compute a numeric jacobian if we're asked to verify the - // analytic jacobian we got - if (_get_jacobian && jacobian_computed && - _sys.verify_analytic_jacobians != 0.0) - { - DenseMatrix analytic_jacobian(_femcontext.elem_jacobian); - - _femcontext.elem_jacobian.zero(); - // Logging of numerical jacobians is done separately - _sys.numerical_side_jacobian(_femcontext); - - Real analytic_norm = analytic_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(_femcontext.elem_jacobian); - - // The matrix "analytic_jacobian" will now hold the error matrix - analytic_jacobian.add(-1.0, _femcontext.elem_jacobian); - Real error_norm = analytic_jacobian.l1_norm(); - - Real relative_error = error_norm / - std::max(analytic_norm, numerical_norm); - - if (relative_error > _sys.verify_analytic_jacobians) - { - libMesh::err << "Relative error " << relative_error - << " detected in analytic jacobian on element " - << _femcontext.elem->id() - << ", side " - << static_cast(_femcontext.side) << '!' << std::endl; - - unsigned int old_precision = libMesh::out.precision(); - libMesh::out.precision(16); - libMesh::out << "J_analytic " << _femcontext.elem->id() << " = " - << _femcontext.elem_jacobian << std::endl; - analytic_jacobian.add(1.0, _femcontext.elem_jacobian); - libMesh::out << "J_numeric " << _femcontext.elem->id() << " = " - << analytic_jacobian << std::endl; - libMesh::out.precision(old_precision); - - libmesh_error(); - } - // Once we've verified a side, we'll want to add back the - // rest of the accumulated 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 -#ifdef DEBUG - if (_get_jacobian && jacobian_computed && - _sys.verify_analytic_jacobians == 0.0) - { - _femcontext.elem_jacobian += old_jacobian; - } -#endif // ifdef DEBUG - } - -#ifdef LIBMESH_ENABLE_CONSTRAINTS - // We turn off the asymmetric constraint application; - // enforce_constraints_exactly() should be called in the solver - if (_get_residual && _get_jacobian) - _sys.get_dof_map().constrain_element_matrix_and_vector - (_femcontext.elem_jacobian, _femcontext.elem_residual, - _femcontext.dof_indices, false); - else if (_get_residual) - _sys.get_dof_map().constrain_element_vector - (_femcontext.elem_residual, _femcontext.dof_indices, false); - else if (_get_jacobian) - _sys.get_dof_map().constrain_element_matrix - (_femcontext.elem_jacobian, _femcontext.dof_indices, false); -#endif // #ifdef LIBMESH_ENABLE_CONSTRAINTS - - if (_get_jacobian && _sys.print_element_jacobians) - { - unsigned int old_precision = libMesh::out.precision(); - libMesh::out.precision(16); - libMesh::out << "J_elem " << _femcontext.elem->id() << " = " - << _femcontext.elem_jacobian << std::endl; - libMesh::out.precision(old_precision); - } - - { // A lock is necessary around access to the global system - femsystem_mutex::scoped_lock lock(assembly_mutex); - - if (_get_jacobian) - _sys.matrix->add_matrix (_femcontext.elem_jacobian, - _femcontext.dof_indices); - if (_get_residual) - _sys.rhs->add_vector (_femcontext.elem_residual, - _femcontext.dof_indices); - } // Scope for assembly mutex - - } - } - - private: - - FEMSystem& _sys; - - const bool _get_residual, _get_jacobian; - }; - - class PostprocessContributions - { - public: - /** - * constructor to set context - */ - explicit - PostprocessContributions(FEMSystem &sys) : _sys(sys) {} - - /** - * operator() for use with Threads::parallel_for(). - */ - void operator()(const ConstElemRange &range) const - { - AutoPtr con = _sys.build_context(); - FEMContext &_femcontext = libmesh_cast_ref(*con); - _sys.init_context(_femcontext); - - for (ConstElemRange::const_iterator elem_it = range.begin(); - elem_it != range.end(); ++elem_it) - { - Elem *el = const_cast(*elem_it); - _femcontext.pre_fe_reinit(_sys, el); - - // Optionally initialize all the interior FE objects on elem. - if (_sys.fe_reinit_during_postprocess) - _femcontext.elem_fe_reinit(); - - _sys.element_postprocess(_femcontext); - - for (_femcontext.side = 0; - _femcontext.side != _femcontext.elem->n_sides(); - ++_femcontext.side) - { - // Don't compute on non-boundary sides unless requested - if (!_sys.postprocess_sides || - (!_sys.compute_internal_sides && - _femcontext.elem->neighbor(_femcontext.side) != NULL)) - continue; - - // Optionally initialize all the FE objects on this side. - if (_sys.fe_reinit_during_postprocess) - _femcontext.side_fe_reinit(); - - _sys.side_postprocess(_femcontext); - } - } - } - - private: - - FEMSystem& _sys; - }; - - class QoIContributions - { - public: - /** - * constructor to set context - */ - explicit - QoIContributions(FEMSystem &sys, DifferentiableQoI &diff_qoi, - const QoISet &qoi_indices) : - qoi(sys.qoi.size(), 0.), _sys(sys), _diff_qoi(diff_qoi),_qoi_indices(qoi_indices) {} - - /** - * splitting constructor - */ - QoIContributions(const QoIContributions &other, - Threads::split) : - qoi(other._sys.qoi.size(), 0.), _sys(other._sys), _diff_qoi(other._diff_qoi) {} - - /** - * operator() for use with Threads::parallel_reduce(). - */ - void operator()(const ConstElemRange &range) - { - AutoPtr con = _sys.build_context(); - FEMContext &_femcontext = libmesh_cast_ref(*con); - _diff_qoi.init_context(_femcontext); - - for (ConstElemRange::const_iterator elem_it = range.begin(); - elem_it != range.end(); ++elem_it) - { - Elem *el = const_cast(*elem_it); - - _femcontext.pre_fe_reinit(_sys, el); - - if (_diff_qoi.assemble_qoi_elements) - { - _femcontext.elem_fe_reinit(); - - _diff_qoi.element_qoi(_femcontext, _qoi_indices); - } - - for (_femcontext.side = 0; - _femcontext.side != _femcontext.elem->n_sides(); - ++_femcontext.side) - { - // Don't compute on non-boundary sides unless requested - if (!_diff_qoi.assemble_qoi_sides || - (!_diff_qoi.assemble_qoi_internal_sides && - _femcontext.elem->neighbor(_femcontext.side) != NULL)) - continue; - - _femcontext.side_fe_reinit(); - - _diff_qoi.side_qoi(_femcontext, _qoi_indices); - } - } - - this->_diff_qoi.thread_join( this->qoi, _femcontext.elem_qoi, _qoi_indices ); - } - - void join (const QoIContributions& other) - { - libmesh_assert_equal_to (this->qoi.size(), other.qoi.size()); - this->_diff_qoi.thread_join( this->qoi, other.qoi, _qoi_indices ); - } - - std::vector qoi; - - private: - - FEMSystem& _sys; - DifferentiableQoI& _diff_qoi; - - const QoISet _qoi_indices; - }; - - class QoIDerivativeContributions - { - public: - /** - * constructor to set context - */ - QoIDerivativeContributions(FEMSystem &sys, const QoISet& qoi_indices, - DifferentiableQoI &qoi ) : - _sys(sys), _qoi_indices(qoi_indices), _qoi(qoi) {} - - /** - * operator() for use with Threads::parallel_for(). - */ - void operator()(const ConstElemRange &range) const - { - AutoPtr con = _sys.build_context(); - FEMContext &_femcontext = libmesh_cast_ref(*con); - _qoi.init_context(_femcontext); - - for (ConstElemRange::const_iterator elem_it = range.begin(); - elem_it != range.end(); ++elem_it) - { - Elem *el = const_cast(*elem_it); - - _femcontext.pre_fe_reinit(_sys, el); - - if (_qoi.assemble_qoi_elements) - { - _femcontext.elem_fe_reinit(); - - _qoi.element_qoi_derivative(_femcontext, _qoi_indices); - } - - for (_femcontext.side = 0; - _femcontext.side != _femcontext.elem->n_sides(); - ++_femcontext.side) - { - // Don't compute on non-boundary sides unless requested - if (!_qoi.assemble_qoi_sides || - (!_qoi.assemble_qoi_internal_sides && - _femcontext.elem->neighbor(_femcontext.side) != NULL)) - continue; - - _femcontext.side_fe_reinit(); - - _qoi.side_qoi_derivative(_femcontext, _qoi_indices); - } - - // We need some unmodified indices to use for constraining - // multiple vector - // FIXME - there should be a DofMap::constrain_element_vectors - // to do this more efficiently - std::vector original_dofs = _femcontext.dof_indices; - - for (unsigned int i=0; i != _sys.qoi.size(); ++i) - if (_qoi_indices.has_index(i)) - { - _femcontext.dof_indices = original_dofs; - _sys.get_dof_map().constrain_element_vector - (_femcontext.elem_qoi_derivative[i], _femcontext.dof_indices, false); - - _sys.get_adjoint_rhs(i).add_vector - (_femcontext.elem_qoi_derivative[i], _femcontext.dof_indices); - } - } - } - - private: - - FEMSystem& _sys; - const QoISet& _qoi_indices; - DifferentiableQoI& _qoi; - }; - - -} - - namespace libMesh { - - - - -FEMSystem::FEMSystem (EquationSystems& es, - const std::string& name, - const unsigned int number) - : Parent(es, name, number), - fe_reinit_during_postprocess(true), - numerical_jacobian_h(TOLERANCE), - verify_analytic_jacobians(0.0) +bool FEMPhysics::eulerian_residual (bool request_jacobian, + DiffContext &/*c*/) { -} - - -FEMSystem::~FEMSystem () -{ - this->clear(); -} - - - -void FEMSystem::clear() -{ - Parent::clear(); -} - - - -void FEMSystem::init_data () -{ - // First initialize LinearImplicitSystem data - Parent::init_data(); -} - - -void FEMSystem::assembly (bool get_residual, bool get_jacobian) -{ - libmesh_assert(get_residual || get_jacobian); - std::string log_name; - if (get_residual && get_jacobian) - log_name = "assembly()"; - else if (get_residual) - log_name = "assembly(get_residual)"; - else - log_name = "assembly(get_jacobian)"; - - START_LOG(log_name, "FEMSystem"); - - const MeshBase& mesh = this->get_mesh(); - -// this->get_vector("_nonlinear_solution").localize -// (*current_local_nonlinear_solution, -// dof_map.get_send_list()); - this->update(); - - if (print_solution_norms) - { -// this->get_vector("_nonlinear_solution").close(); - this->solution->close(); - - unsigned int old_precision = libMesh::out.precision(); - libMesh::out.precision(16); - libMesh::out << "|U| = " -// << this->get_vector("_nonlinear_solution").l1_norm() - << this->solution->l1_norm() - << std::endl; - libMesh::out.precision(old_precision); - } - if (print_solutions) - { - unsigned int old_precision = libMesh::out.precision(); - libMesh::out.precision(16); -// libMesh::out << "U = [" << this->get_vector("_nonlinear_solution") - libMesh::out << "U = [" << *(this->solution) - << "];" << std::endl; - libMesh::out.precision(old_precision); - } - - // Is this definitely necessary? [RHS] - // Yes. [RHS 2012] - if (get_jacobian) - matrix->zero(); - if (get_residual) - rhs->zero(); - - // Stupid C++ lets you set *Real* verify_analytic_jacobians = true! - if (verify_analytic_jacobians > 0.5) - { - libMesh::err << "WARNING! verify_analytic_jacobians was set " - << "to absurdly large value of " - << verify_analytic_jacobians << std::endl; - libMesh::err << "Resetting to 1e-6!" << std::endl; - verify_analytic_jacobians = 1e-6; - } - - // In time-dependent problems, the nonlinear function we're trying - // to solve at each timestep may depend on the particular solver - // we're using - libmesh_assert(time_solver.get()); - - // Build the residual and jacobian contributions on every active - // mesh element on this processor - Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(), - mesh.active_local_elements_end()), - AssemblyContributions(*this, get_residual, get_jacobian)); - - - if (get_residual && (print_residual_norms || print_residuals)) - this->rhs->close(); - if (get_residual && print_residual_norms) - { - unsigned int old_precision = libMesh::out.precision(); - libMesh::out.precision(16); - libMesh::out << "|F| = " << this->rhs->l1_norm() << std::endl; - libMesh::out.precision(old_precision); - } - if (get_residual && print_residuals) - { - unsigned int old_precision = libMesh::out.precision(); - libMesh::out.precision(16); - libMesh::out << "F = [" << *(this->rhs) << "];" << std::endl; - libMesh::out.precision(old_precision); - } - - if (get_jacobian && (print_jacobian_norms || print_jacobians)) - this->matrix->close(); - if (get_jacobian && print_jacobian_norms) - { - unsigned int old_precision = libMesh::out.precision(); - libMesh::out.precision(16); - libMesh::out << "|J| = " << this->matrix->l1_norm() << std::endl; - libMesh::out.precision(old_precision); - } - if (get_jacobian && print_jacobians) - { - unsigned int old_precision = libMesh::out.precision(); - libMesh::out.precision(16); - libMesh::out << "J = [" << *(this->matrix) << "];" << std::endl; - libMesh::out.precision(old_precision); - } - STOP_LOG(log_name, "FEMSystem"); -} - - - -void FEMSystem::solve() -{ - Parent::solve(); - - // On a moving mesh we want the mesh to reflect the new solution - this->mesh_position_set(); -} - - - -void FEMSystem::mesh_position_set() -{ - // If we don't need to move the mesh, we're done - if (_mesh_sys != this) - return; - - const MeshBase& mesh = this->get_mesh(); - - // This code won't work on a parallelized mesh yet - - // it won't get ancestor elements right. - libmesh_assert(mesh.is_serial()); - - // In fact it won't work in parallel at all yet - we don't have - // current_solution() data for non-ghost elements. That's going to - // be put on hold until we're debugged in serial and until I figure - // out how to better abstract the "ask other procs for their local - // data, respond to others' queries, work on my query's results" - // pattern that we seem to be using a lot. - libmesh_assert_equal_to (libMesh::n_processors(), 1); - - AutoPtr con = this->build_context(); - FEMContext &_femcontext = libmesh_cast_ref(*con); - this->init_context(_femcontext); - - // Move every mesh element we can - MeshBase::const_element_iterator el = - mesh.active_elements_begin(); - const MeshBase::const_element_iterator end_el = - mesh.active_elements_end(); - - for ( ; el != end_el; ++el) - { - // We need the algebraic data - _femcontext.pre_fe_reinit(*this, *el); - // And when asserts are on, we also need the FE so - // we can assert that the mesh data is of the right type. -#ifndef NDEBUG - _femcontext.elem_fe_reinit(); -#endif - - // This code won't handle moving subactive elements - libmesh_assert(!_femcontext.elem->has_children()); - - _femcontext.elem_position_set(1.); - } -} - - - -void FEMSystem::postprocess () -{ - START_LOG("postprocess()", "FEMSystem"); - - const MeshBase& mesh = this->get_mesh(); - - this->update(); - - // Loop over every active mesh element on this processor - Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(), - mesh.active_local_elements_end()), - PostprocessContributions(*this)); - - STOP_LOG("postprocess()", "FEMSystem"); -} - - - -void FEMSystem::assemble_qoi (const QoISet &qoi_indices) -{ - START_LOG("assemble_qoi()", "FEMSystem"); - - const MeshBase& mesh = this->get_mesh(); - - this->update(); - - // the quantity of interest is assumed to be a sum of element and - // side terms - for (unsigned int i=0; i != qoi.size(); ++i) - if (qoi_indices.has_index(i)) - qoi[i] = 0; - - // Create a non-temporary qoi_contributions object, so we can query - // its results after the reduction - QoIContributions qoi_contributions(*this, *(this->diff_qoi), qoi_indices); - - // Loop over every active mesh element on this processor - Threads::parallel_reduce(elem_range.reset(mesh.active_local_elements_begin(), - mesh.active_local_elements_end()), - qoi_contributions); - - this->diff_qoi->parallel_op( this->qoi, qoi_contributions.qoi, qoi_indices ); - - STOP_LOG("assemble_qoi()", "FEMSystem"); -} - - - -void FEMSystem::assemble_qoi_derivative (const QoISet& qoi_indices) -{ - START_LOG("assemble_qoi_derivative()", "FEMSystem"); - - const MeshBase& mesh = this->get_mesh(); - - this->update(); - - // The quantity of interest derivative assembly accumulates on - // initially zero vectors - for (unsigned int i=0; i != qoi.size(); ++i) - if (qoi_indices.has_index(i)) - this->add_adjoint_rhs(i).zero(); - - // Loop over every active mesh element on this processor - Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(), - mesh.active_local_elements_end()), - QoIDerivativeContributions(*this, qoi_indices, - *(this->diff_qoi))); - - STOP_LOG("assemble_qoi_derivative()", "FEMSystem"); -} - - - -void FEMSystem::numerical_jacobian (TimeSolverResPtr res, - FEMContext &context) -{ - // Logging is done by numerical_elem_jacobian - // or numerical_side_jacobian - - DenseVector original_residual(context.elem_residual); - DenseVector backwards_residual(context.elem_residual); - DenseMatrix numerical_jacobian(context.elem_jacobian); -#ifdef DEBUG - DenseMatrix old_jacobian(context.elem_jacobian); -#endif - - Real numerical_point_h = 0.; - if (_mesh_sys == this) - numerical_point_h = numerical_jacobian_h * context.elem->hmin(); - - for (unsigned int j = 0; j != context.dof_indices.size(); ++j) - { - // Take the "minus" side of a central differenced first derivative - 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 - Real *coord = NULL; - if (_mesh_sys == this) - { - if (_mesh_x_var != libMesh::invalid_uint) - for (unsigned int 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 = &(const_cast(context.elem)->point(k)(0)); - if (_mesh_y_var != libMesh::invalid_uint) - for (unsigned int k = 0; - k != context.dof_indices_var[_mesh_y_var].size(); ++k) - if (context.dof_indices_var[_mesh_y_var][k] == - context.dof_indices[j]) - coord = &(const_cast(context.elem)->point(k)(1)); - if (_mesh_z_var != libMesh::invalid_uint) - for (unsigned int k = 0; - k != context.dof_indices_var[_mesh_z_var].size(); ++k) - if (context.dof_indices_var[_mesh_z_var][k] == - context.dof_indices[j]) - coord = &(const_cast(context.elem)->point(k)(2)); - } - if (coord) - { - // We have enough information to scale the perturbations - // here appropriately - context.elem_solution(j) = original_solution - numerical_point_h; - *coord = libmesh_real(context.elem_solution(j)); - } - - context.elem_residual.zero(); - ((*time_solver).*(res))(false, context); -#ifdef DEBUG - libmesh_assert_equal_to (old_jacobian, context.elem_jacobian); -#endif - backwards_residual = context.elem_residual; - - // Take the "plus" side of a central differenced first derivative - context.elem_solution(j) = original_solution + numerical_jacobian_h; - if (coord) - { - context.elem_solution(j) = original_solution + numerical_point_h; - *coord = libmesh_real(context.elem_solution(j)); - } - context.elem_residual.zero(); - ((*time_solver).*(res))(false, context); -#ifdef DEBUG - libmesh_assert_equal_to (old_jacobian, context.elem_jacobian); -#endif - - context.elem_solution(j) = original_solution; - if (coord) - { - *coord = libmesh_real(context.elem_solution(j)); - for (unsigned int i = 0; i != context.dof_indices.size(); ++i) - { - numerical_jacobian(i,j) = - (context.elem_residual(i) - backwards_residual(i)) / - 2. / numerical_point_h; - } - } - else - { - for (unsigned int i = 0; i != context.dof_indices.size(); ++i) - { - numerical_jacobian(i,j) = - (context.elem_residual(i) - backwards_residual(i)) / - 2. / numerical_jacobian_h; - } - } - } - - context.elem_residual = original_residual; - context.elem_jacobian = numerical_jacobian; -} - - - -void FEMSystem::numerical_elem_jacobian (FEMContext &context) -{ - START_LOG("numerical_elem_jacobian()", "FEMSystem"); - this->numerical_jacobian(&TimeSolver::element_residual, context); - STOP_LOG("numerical_elem_jacobian()", "FEMSystem"); -} - - - -void FEMSystem::numerical_side_jacobian (FEMContext &context) -{ - START_LOG("numerical_side_jacobian()", "FEMSystem"); - this->numerical_jacobian(&TimeSolver::side_residual, context); - STOP_LOG("numerical_side_jacobian()", "FEMSystem"); -} - - - -AutoPtr FEMSystem::build_context () -{ - AutoPtr ap(new FEMContext(*this)); - - ap->set_deltat_pointer( &deltat ); - - return ap; -} - - - -void FEMSystem::init_context(DiffContext &c) -{ - // Parent::init_context(c); // may be a good idea in derived classes - - // Although we do this in DiffSystem::build_context() and - // FEMSystem::build_context() as well, we do it here just to be - // extra sure that the deltat pointer gets set. Since the - // intended behavior is for classes derived from FEMSystem to - // call Parent::init_context() in their own init_context() - // overloads, we can ensure that those classes get the correct - // deltat pointers even if they have different build_context() - // overloads. - c.set_deltat_pointer ( &deltat ); - - 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); -} - - - -void FEMSystem::mesh_position_get() -{ - // This function makes no sense unless we've already picked out some - // variable(s) to reflect mesh position coordinates - if (!_mesh_sys) - libmesh_error(); - - // We currently assume mesh variables are in our own system - if (_mesh_sys != this) - libmesh_not_implemented(); - - // Loop over every active mesh element on this processor - const MeshBase& mesh = this->get_mesh(); - - MeshBase::const_element_iterator el = - mesh.active_local_elements_begin(); - 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) - { - _femcontext.pre_fe_reinit(*this, *el); - - _femcontext.elem_position_get(); - - if (_mesh_x_var != libMesh::invalid_uint) - this->solution->insert(*_femcontext.elem_subsolutions[_mesh_x_var], - _femcontext.dof_indices_var[_mesh_x_var]); - if (_mesh_y_var != libMesh::invalid_uint) - this->solution->insert(*_femcontext.elem_subsolutions[_mesh_y_var], - _femcontext.dof_indices_var[_mesh_y_var]); - if (_mesh_z_var != libMesh::invalid_uint) - this->solution->insert(*_femcontext.elem_subsolutions[_mesh_z_var], - _femcontext.dof_indices_var[_mesh_z_var]); - } - - this->solution->close(); - - // And make sure the current_local_solution is up to date too - this->update(); -} - - -bool FEMSystem::eulerian_residual (bool request_jacobian, - DiffContext &c) -{ // Only calculate a mesh movement residual if it's necessary if (!_mesh_sys) return request_jacobian; + libmesh_not_implemented(); + +#if 0 FEMContext &context = libmesh_cast_ref(c); // This function only supports fully coupled mesh motion for now @@ -1027,15 +79,15 @@ const std::vector > &psi = context.element_fe_var[mesh_xyz_var]->get_phi(); - for (unsigned int var = 0; var != this->n_vars(); ++var) + for (unsigned int var = 0; var != context.n_vars(); ++var) { // Mesh motion only affects time-evolving variables - if (!_time_evolving[var]) + if (this->is_time_evolving(var)) continue; // The mesh coordinate variables themselves are Lagrangian, // not Eulerian, and no convective term is desired. - if (_mesh_sys == this && + if (/*_mesh_sys == this && */ (var == _mesh_x_var || var == _mesh_y_var || var == _mesh_z_var)) @@ -1139,22 +191,23 @@ } } } +#endif // 0 return request_jacobian; } -bool FEMSystem::mass_residual (bool request_jacobian, - DiffContext &c) +bool FEMPhysics::mass_residual (bool request_jacobian, + DiffContext &c) { FEMContext &context = libmesh_cast_ref(c); unsigned int n_qpoints = (context.get_element_qrule())->n_points(); - for (unsigned int var = 0; var != this->n_vars(); ++var) + for (unsigned int var = 0; var != context.n_vars(); ++var) { - if (!_time_evolving[var]) + if (this->is_time_evolving(var)) continue; const std::vector &JxW = Index: src/systems/system.C =================================================================== --- src/systems/system.C (revision 6309) +++ src/systems/system.C (working copy) @@ -64,10 +64,6 @@ current_local_solution (NumericVector::build()), time (0.), qoi (0), - _mesh_sys (NULL), - _mesh_x_var (libMesh::invalid_uint), - _mesh_y_var (libMesh::invalid_uint), - _mesh_z_var (libMesh::invalid_uint), _init_system_function (NULL), _init_system_object (NULL), _assemble_system_function (NULL), Index: src/systems/fem_system.C =================================================================== --- src/systems/fem_system.C (revision 6309) +++ src/systems/fem_system.C (working copy) @@ -135,7 +135,7 @@ ++_femcontext.side) { // Don't compute on non-boundary sides unless requested - if (!_sys.compute_internal_sides && + if (!_sys.get_physics()->compute_internal_sides && _femcontext.elem->neighbor(_femcontext.side) != NULL) continue; @@ -319,7 +319,7 @@ { // Don't compute on non-boundary sides unless requested if (!_sys.postprocess_sides || - (!_sys.compute_internal_sides && + (!_sys.get_physics()->compute_internal_sides && _femcontext.elem->neighbor(_femcontext.side) != NULL)) continue; @@ -922,7 +922,7 @@ // Make sure we're prepared to do mass integration for (unsigned int var = 0; var != this->n_vars(); ++var) - if (_time_evolving[var]) + if (this->get_physics()->is_time_evolving(var)) { context.element_fe_var[var]->get_JxW(); context.element_fe_var[var]->get_phi(); @@ -931,14 +931,6 @@ -void FEMSystem::time_evolving (unsigned int var) -{ - // Call the parent function - Parent::time_evolving(var); -} - - - void FEMSystem::mesh_position_get() { // This function makes no sense unless we've already picked out some @@ -986,213 +978,4 @@ this->update(); } - -bool FEMSystem::eulerian_residual (bool request_jacobian, - DiffContext &c) -{ - // Only calculate a mesh movement residual if it's necessary - if (!_mesh_sys) - return request_jacobian; - - FEMContext &context = libmesh_cast_ref(c); - - // This function only supports fully coupled mesh motion for now - libmesh_assert_equal_to (_mesh_sys, this); - - unsigned int n_qpoints = (context.get_element_qrule())->n_points(); - - const unsigned int n_x_dofs = (_mesh_x_var == libMesh::invalid_uint) ? - 0 : context.dof_indices_var[_mesh_x_var].size(); - const unsigned int n_y_dofs = (_mesh_y_var == libMesh::invalid_uint) ? - 0 : context.dof_indices_var[_mesh_y_var].size(); - const unsigned int n_z_dofs = (_mesh_z_var == libMesh::invalid_uint) ? - 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 : - (n_z_dofs ? _mesh_z_var : - libMesh::invalid_uint)); - - // If we're our own _mesh_sys, we'd better be in charge of - // at least one coordinate, and we'd better have the same - // FE type for all coordinates we are in charge of - libmesh_assert_not_equal_to (mesh_xyz_var, libMesh::invalid_uint); - 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 = - context.element_fe_var[mesh_xyz_var]->get_phi(); - - for (unsigned int var = 0; var != this->n_vars(); ++var) - { - // Mesh motion only affects time-evolving variables - if (!_time_evolving[var]) - continue; - - // The mesh coordinate variables themselves are Lagrangian, - // not Eulerian, and no convective term is desired. - if (_mesh_sys == this && - (var == _mesh_x_var || - var == _mesh_y_var || - var == _mesh_z_var)) - continue; - - // Some of this code currently relies on the assumption that - // we can pull mesh coordinate data from our own system - if (_mesh_sys != this) - libmesh_not_implemented(); - - // This residual should only be called by unsteady solvers: - // if the mesh is steady, there's no mesh convection term! - UnsteadySolver *unsteady; - if (this->time_solver->is_steady()) - return request_jacobian; - else - unsteady = libmesh_cast_ptr(this->time_solver.get()); - - const std::vector &JxW = - context.element_fe_var[var]->get_JxW(); - - const std::vector > &phi = - context.element_fe_var[var]->get_phi(); - - const std::vector > &dphi = - context.element_fe_var[var]->get_dphi(); - - const unsigned int n_u_dofs = context.dof_indices_var[var].size(); - - DenseSubVector &Fu = *context.elem_subresiduals[var]; - DenseSubMatrix &Kuu = *context.elem_subjacobians[var][var]; - - DenseSubMatrix *Kux = n_x_dofs ? - context.elem_subjacobians[var][_mesh_x_var] : NULL; - DenseSubMatrix *Kuy = n_y_dofs ? - context.elem_subjacobians[var][_mesh_y_var] : NULL; - DenseSubMatrix *Kuz = n_z_dofs ? - context.elem_subjacobians[var][_mesh_z_var] : NULL; - - std::vector delta_x(n_x_dofs, 0.); - std::vector delta_y(n_y_dofs, 0.); - std::vector delta_z(n_z_dofs, 0.); - - for (unsigned int i = 0; i != n_x_dofs; ++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 = 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 = 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 = context.interior_gradient(var, qp); - RealGradient convection(0.); - - for (unsigned int i = 0; i != n_x_dofs; ++i) - convection(0) += delta_x[i] * psi[i][qp]; - for (unsigned int i = 0; i != n_y_dofs; ++i) - convection(1) += delta_y[i] * psi[i][qp]; - for (unsigned int i = 0; i != n_z_dofs; ++i) - convection(2) += delta_z[i] * psi[i][qp]; - - for (unsigned int i = 0; i != n_u_dofs; ++i) - { - Number JxWxPhiI = JxW[qp] * phi[i][qp]; - Fu(i) += (convection * grad_u) * JxWxPhiI; - if (request_jacobian) - { - Number JxWxPhiI = JxW[qp] * phi[i][qp]; - for (unsigned int j = 0; j != n_u_dofs; ++j) - Kuu(i,j) += JxWxPhiI * (convection * dphi[j][qp]); - - Number JxWxPhiIoverDT = JxWxPhiI/this->deltat; - - Number JxWxPhiIxDUDXoverDT = JxWxPhiIoverDT * grad_u(0); - for (unsigned int j = 0; j != n_x_dofs; ++j) - (*Kux)(i,j) += JxWxPhiIxDUDXoverDT * psi[j][qp]; - - Number JxWxPhiIxDUDYoverDT = JxWxPhiIoverDT * grad_u(1); - for (unsigned int j = 0; j != n_y_dofs; ++j) - (*Kuy)(i,j) += JxWxPhiIxDUDYoverDT * psi[j][qp]; - - Number JxWxPhiIxDUDZoverDT = JxWxPhiIoverDT * grad_u(2); - for (unsigned int j = 0; j != n_z_dofs; ++j) - (*Kuz)(i,j) += JxWxPhiIxDUDZoverDT * psi[j][qp]; - } - } - } - } - - return request_jacobian; -} - - - -bool FEMSystem::mass_residual (bool request_jacobian, - DiffContext &c) -{ - FEMContext &context = libmesh_cast_ref(c); - - unsigned int n_qpoints = (context.get_element_qrule())->n_points(); - - for (unsigned int var = 0; var != this->n_vars(); ++var) - { - if (!_time_evolving[var]) - continue; - - const std::vector &JxW = - context.element_fe_var[var]->get_JxW(); - - const std::vector > &phi = - context.element_fe_var[var]->get_phi(); - - const unsigned int n_dofs = context.dof_indices_var[var].size(); - - DenseSubVector &Fu = *context.elem_subresiduals[var]; - DenseSubMatrix &Kuu = *context.elem_subjacobians[var][var]; - - for (unsigned int qp = 0; qp != n_qpoints; ++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 && context.elem_solution_derivative) - { - libmesh_assert_equal_to (context.elem_solution_derivative, 1.0); - - Number JxWxPhiI = JxW[qp] * phi[i][qp]; - Kuu(i,i) += JxWxPhiI * phi[i][qp]; - for (unsigned int j = i+1; j < n_dofs; ++j) - { - Number Kij = JxWxPhiI * phi[j][qp]; - Kuu(i,j) += Kij; - Kuu(j,i) += Kij; - } - } - } - } - } - - return request_jacobian; -} - } // namespace libMesh Index: src/systems/fem_context.C =================================================================== --- src/systems/fem_context.C (revision 6309) +++ src/systems/fem_context.C (working copy) @@ -37,10 +37,10 @@ : DiffContext(sys), element_qrule(NULL), side_qrule(NULL), edge_qrule(NULL), - _mesh_sys(sys.get_mesh_system()), - _mesh_x_var(sys.get_mesh_x_var()), - _mesh_y_var(sys.get_mesh_y_var()), - _mesh_z_var(sys.get_mesh_z_var()), + _mesh_sys(NULL), + _mesh_x_var(0), + _mesh_y_var(0), + _mesh_z_var(0), elem(NULL), side(0), edge(0), dim(sys.get_mesh().mesh_dimension()), _boundary_info(sys.get_mesh().boundary_info.get()) Index: src/systems/diff_system.C =================================================================== --- src/systems/diff_system.C (revision 6309) +++ src/systems/diff_system.C (working copy) @@ -39,7 +39,7 @@ print_jacobian_norms(false), print_jacobians(false), print_element_jacobians(false), - diff_physics(this), + _diff_physics(this), diff_qoi(this) { } @@ -55,7 +55,15 @@ void DifferentiableSystem::clear () { - this->clear_physics(); + // If we had an attached Physics object, delete it. + if (this->_diff_physics != this) + { + delete this->_diff_physics; + this->_diff_physics = this; + } + // If we had no attached Physics object, clear our own Physics data + else + this->clear_physics(); // If we had an attached QoI object, delete it. if (this->diff_qoi != this) @@ -86,8 +94,10 @@ void DifferentiableSystem::init_data () { - // Do any initialization our physics needs - this->init_physics(*this); + // If it isn't a separate initialized-upon-attachment object, do any + // initialization our physics needs. + if (this->_diff_physics == this) + this->init_physics(*this); // Do any initialization our solvers need libmesh_assert(time_solver.get());