cl_TimestepMatrices.hpp

Go to the documentation of this file.

/*
 * BELFEM -- The Berkeley Lab Finite Element Framework
 * Copyright (c) 2026, The Regents of the University of California, through
 * Lawrence Berkeley National Laboratory (subject to receipt of any required
 * approvals from the U.S. Dept. of Energy).  All rights reserved.
 *
 * Developers: Christian Messe, Gregory Giard
 *
 * See the top-level LICENSE file for the complete license and disclaimer.
 */
 
#ifndef BELFEM_CL_FEM_TIMESTEPMATRICES_HPP
#define BELFEM_CL_FEM_TIMESTEPMATRICES_HPP
 
#include "typedefs.hpp"
#include "cl_Vector.hpp"
#include "cl_Matrix.hpp"
#include "cl_Bitset.hpp"
 
namespace belfem
{
    namespace fem
    {
        enum class MatrixFlag : index_t
        {
            M            = 0,   
            D            = 1,   
            K            = 2,   
            F            = 3,   
            dMdX_times_x = 4,   
            dMdX_times_h = 5,   
            dKdX_times_x = 6,   
            dFdX         = 7    
        };

        class TimestepMatrices
        {
            //--------------------------------------------------------------
            // Physics matrices (populated by IWG) - all DENSE, n_e × n_e
            //--------------------------------------------------------------
 
            Matrix< real > mM;    
            Matrix< real > mD;    
            Matrix< real > mK;    
            Vector< real > mF;    
 
            // Contracted derivatives: (dX_ik/dx_j) * v_k for some vector v
            // These are dense n_e × n_e matrices, computed by IWG
            Matrix< real > mdMdX_times_x;   
            Matrix< real > mdMdX_times_h;   
            Matrix< real > mdKdX_times_x;   
            Matrix< real > mdFdX;           
 
            //--------------------------------------------------------------
            // Assembled element system (computed by time-stepping scheme)
            //--------------------------------------------------------------
 
            Matrix< real > mJ;    
            Matrix< real > mdJdx;    
 
            //--------------------------------------------------------------
            // Flags
            //--------------------------------------------------------------
 
            Bitset< 8 > mFlags;
 
        public:
 
            //--------------------------------------------------------------
            // Constructor / Destructor
            //--------------------------------------------------------------
 
            TimestepMatrices() = default;
            ~TimestepMatrices() = default;

            void
            initialize( const index_t aNumDofs );
 
            void
            reset();
 
            //--------------------------------------------------------------
            // Flag management
            //--------------------------------------------------------------

            void
            reset_flags()
            {
                mFlags.reset();
            }
 
            void
            set_flag( const MatrixFlag aFlag )
            {
                mFlags.set( static_cast< index_t >( aFlag ) );
            }
 
            bool
            has_flag( const MatrixFlag aFlag ) const
            {
                return mFlags.test( static_cast< index_t >( aFlag ) );
            }

            const Bitset< 8 > &
            flags() const
            {
                return mFlags;
            }
 
            //--------------------------------------------------------------
            // Const accessors
            //--------------------------------------------------------------
 
            const Matrix< real > & M() const { return mM; }
            const Matrix< real > & D() const { return mD; }
            const Matrix< real > & K() const { return mK; }
            const Vector< real > & f() const { return mF; }
 
            const Matrix< real > & dMdx_times_x() const { return mdMdX_times_x; }
            const Matrix< real > & dMdx_times_h() const { return mdMdX_times_h; }
            const Matrix< real > & dKdx_times_x() const { return mdKdX_times_x; }
            const Matrix< real > & dfdx()         const { return mdFdX; }
 
            const Matrix< real > & J() const { return mJ; }
            const Matrix< real > & dJdx() const { return mdJdx; }
 
            //--------------------------------------------------------------
            // Non-const accessors (for IWG to fill)
            //--------------------------------------------------------------
 
            Matrix< real > & M() { return mM; }
            Matrix< real > & D() { return mD; }
            Matrix< real > & K() { return mK; }
            Vector< real > & f() { return mF; }
 
            Matrix< real > & dMdx_times_x() { return mdMdX_times_x; }
            Matrix< real > & dMdx_times_h() { return mdMdX_times_h; }
            Matrix< real > & dKdx_times_x() { return mdKdX_times_x; }
            Matrix< real > & dfdx()         { return mdFdX; }
 
            //--------------------------------------------------------------
            // Assembly methods
            //--------------------------------------------------------------
 
            void
            assemble_J( const real adt );

            void
            assemble_dJdx( const real adt, const real aAlpha = 1.0 );
 
        };
 
    } // namespace fem
} // namespace belfem
 
#endif // BELFEM_CL_FEM_TIMESTEPMATRICES_HPP