<commit_before>///////////////////////////////////////////////////////////////////////////////
//
// File LinearElasticSystem.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: LinearElasticSystem solve routines 
//
///////////////////////////////////////////////////////////////////////////////

#include <LocalRegions/MatrixKey.h>
#include <MultiRegions/ContField2D.h>
#include <MultiRegions/ContField3D.h>
#include <MultiRegions/GlobalLinSysDirectStaticCond.h>
#include <MultiRegions/GlobalLinSysIterativeStaticCond.h>
#include <MultiRegions/Preconditioner.h>
#include <LinearElasticSolver/EquationSystems/LinearElasticSystem.h>

namespace Nektar
{
    string LinearElasticSystem::className = GetEquationSystemFactory().
        RegisterCreatorFunction("LinearElasticSystem",
                                LinearElasticSystem::create);

    LinearElasticSystem::LinearElasticSystem(
            const LibUtilities::SessionReaderSharedPtr& pSession)
        : EquationSystem(pSession)
    {
    }

    /**
     * @brief Set up the linear elasticity system.
     *
     * This routine loads the E and nu variables from the session file, creates
     * a coupled assembly map and creates containers for the statically
     * condensed block matrix system.
     */
    void LinearElasticSystem::v_InitObject()
    {
        EquationSystem::v_InitObject();

        const int nVel = m_fields[0]->GetCoordim(0);
        int n;

        ASSERTL0(nVel > 1, "Linear elastic solver not set up for"
                           " this dimension (only 2D/3D supported).");

        // Make sure that we have Young's modulus and Poisson ratio set.
        m_session->LoadParameter("E", m_E, 1.0);
        m_session->LoadParameter("nu", m_nu, 0.25);
        m_session->LoadParameter("Beta", m_beta, 1.0e-04);

        // Create a coupled assembly map which allows us to tie u, v and w
        // fields together.
        if (nVel == 2)
        {
            MultiRegions::ContField2DSharedPtr u = boost::dynamic_pointer_cast<
                MultiRegions::ContField2D>(m_fields[0]);
            m_assemblyMap = MemoryManager<CoupledAssemblyMap>
                ::AllocateSharedPtr(m_session,
                                    m_graph,
                                    u->GetLocalToGlobalMap(),
                                    m_boundaryConditions,
                                    m_fields);
        }
        
        if (nVel == 3)
        {
            MultiRegions::ContField3DSharedPtr u = boost::dynamic_pointer_cast<
                MultiRegions::ContField3D>(m_fields[0]);
            m_assemblyMap = MemoryManager<CoupledAssemblyMap>
                ::AllocateSharedPtr(m_session,
                                    m_graph,
                                    u->GetLocalToGlobalMap(),
                                    m_boundaryConditions,
                                    m_fields);
        }
        
        // Figure out size of our new matrix systems by looping over all
        // expansions and multiply number of coefficients by velocity
        // components.
        const int nEl = m_fields[0]->GetExpSize();
        LocalRegions::ExpansionSharedPtr exp;

        Array<OneD, unsigned int> sizeBnd(nEl);
        Array<OneD, unsigned int> sizeInt(nEl);

        // Allocate storage for boundary and interior map caches.
        m_bmap = Array<OneD, Array<OneD, unsigned int> >(nEl);
        m_imap = Array<OneD, Array<OneD, unsigned int> >(nEl);

        // Cache interior and boundary maps to avoid recalculating
        // constantly. Should really be handled by manager in StdRegions but not
        // really working yet.
        for (n = 0; n < nEl; ++n)
        {
            exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
            sizeBnd[n] = nVel * exp->NumBndryCoeffs();
            sizeInt[n] = nVel * exp->GetNcoeffs() - sizeBnd[n];
            exp->GetBoundaryMap(m_bmap[n]);
            exp->GetInteriorMap(m_imap[n]);
        }

        // Create block matrix storage for the statically condensed system.
        MatrixStorage blkmatStorage = eDIAGONAL;
        m_schurCompl = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            sizeBnd, sizeBnd, blkmatStorage);
        m_BinvD      = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            sizeBnd, sizeInt, blkmatStorage);
        m_C          = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            sizeInt, sizeBnd, blkmatStorage);
        m_Dinv       = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            sizeInt, sizeInt, blkmatStorage);
    }

    /**
     * @brief Build matrix system for linear elasticity equations.
     *
     * This routine constructs the matrix discretisation arising from the weak
     * formulation of the linear elasticity equations. The resulting matrix
     * system is then passed to LinearElasticSystem::SetStaticCondBlock in order
     * to construct the statically condensed system.
     *
     * All of the matrices involved in the construction of the divergence of the
     * stress tensor are Laplacian-like. We use the variable coefficient
     * functionality within the library to multiply the Laplacian by the
     * appropriate constants for all diagonal terms, and off-diagonal terms use
     * LinearElasticSystem::BuildLaplacianIJMatrix to construct matrices
     * containing
     *
     * \f[ \int \partial_{x_i} \phi_k \partial_{x_j} \phi_l \f]
     *
     * where mixed derivative terms are present. Symmetry (in terms of k,l) is
     * exploited to avoid constructing this matrix repeatedly.
     *
     * @todo Make static condensation optional and construct full system
     *       instead.
     */
    void LinearElasticSystem::BuildMatrixSystem()
    {
        const int nEl = m_fields[0]->GetExpSize();
        const int nVel = m_fields[0]->GetCoordim(0);
        
        LocalRegions::ExpansionSharedPtr exp;
        int n;

        // Factors map for matrix keys.
        StdRegions::ConstFactorMap factors;

        // Calculate various constants

        // mu
        NekDouble a = m_E * (1.0 - m_nu) / (1.0 + m_nu) / (1.0 - 2.0*m_nu);
        // lambda
        NekDouble b = m_E * 0.5 / (1.0 + m_nu);
        // lambda + mu?
        NekDouble c = m_E * m_nu / (1.0 + m_nu) / (1.0 - 2.0*m_nu);

        bool verbose = m_session->DefinesCmdLineArgument("verbose");
        bool root    = m_session->GetComm()->GetRank() == 0;

        // Loop over each element and construct matrices.
        if (nVel == 2)
        {
            for (n = 0; n < nEl; ++n)
            {
                exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
                const int nPhys = exp->GetTotPoints();
                Array<OneD, NekDouble> aArr(nPhys, a);
                Array<OneD, NekDouble> bArr(nPhys, b);
                
                StdRegions::VarCoeffMap varcoeffA, varcoeffD;
                varcoeffA[StdRegions::eVarCoeffD00] = aArr;
                varcoeffA[StdRegions::eVarCoeffD11] = bArr;
                varcoeffD[StdRegions::eVarCoeffD00] = bArr;
                varcoeffD[StdRegions::eVarCoeffD11] = aArr;

                LocalRegions::MatrixKey matkeyA(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffA);
                LocalRegions::MatrixKey matkeyD(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffD);

                /*
                 * mat holds the linear operator [ A B ] acting on [ u ].
                 *                               [ C D ]           [ v ]
                 */
                Array<TwoD, DNekMatSharedPtr> mat(2,2);
                mat[0][0] = exp->GenMatrix(matkeyA);
                mat[0][1] = BuildLaplacianIJMatrix(1, 0, c, exp);
                mat[1][0] = mat[0][1];
                mat[1][1] = exp->GenMatrix(matkeyD);

                // Set up the statically condensed block for this element.
                SetStaticCondBlock(n, exp, mat);

                if (verbose && root)
                {
                    cout << "\rBuilding matrix system: "
                         << (int)(100.0 * n / nEl) << "%" << flush;
                }
            }
        }
        else if (nVel == 3)
        {
            for (n = 0; n < nEl; ++n)
            {
                exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
                const int nPhys = exp->GetTotPoints();
                Array<OneD, NekDouble> aArr(nPhys, a);
                Array<OneD, NekDouble> bArr(nPhys, b);

                StdRegions::VarCoeffMap varcoeffA, varcoeffE, varcoeffI;
                varcoeffA[StdRegions::eVarCoeffD00] = aArr;
                varcoeffA[StdRegions::eVarCoeffD11] = bArr;
                varcoeffA[StdRegions::eVarCoeffD22] = bArr;
                varcoeffE[StdRegions::eVarCoeffD00] = bArr;
                varcoeffE[StdRegions::eVarCoeffD11] = aArr;
                varcoeffE[StdRegions::eVarCoeffD22] = bArr;
                varcoeffI[StdRegions::eVarCoeffD00] = bArr;
                varcoeffI[StdRegions::eVarCoeffD11] = bArr;
                varcoeffI[StdRegions::eVarCoeffD22] = aArr;

                LocalRegions::MatrixKey matkeyA(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffA);
                LocalRegions::MatrixKey matkeyE(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffE);
                LocalRegions::MatrixKey matkeyI(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffI);

                /*
                 * mat holds the linear operator [ A B C ] acting on [ u ].
                 *                               [ D E F ]           [ v ]
                 *                               [ G H I ]           [ w ]
                 */
                Array<TwoD, DNekMatSharedPtr> mat(3,3);
                mat[0][0] = exp->GenMatrix(matkeyA);
                mat[0][1] = BuildLaplacianIJMatrix(1, 0, c, exp);
                mat[0][2] = BuildLaplacianIJMatrix(2, 0, c, exp);
                
                mat[1][0] = mat[0][1];
                mat[1][1] = exp->GenMatrix(matkeyE);
                mat[1][2] = BuildLaplacianIJMatrix(2, 1, c, exp);
                
                mat[2][0] = mat[0][2];
                mat[2][1] = mat[1][2];
                mat[2][2] = exp->GenMatrix(matkeyI);
                
                // Set up the statically condensed block for this element.
                SetStaticCondBlock(n, exp, mat);

                if (verbose && root)
                {
                    cout << "\rBuilding matrix system: "
                         << (int)(100.0 * n / nEl) << "%" << flush;
                }
            }
        }

        if (verbose && root)
        {
            cout << "\rBuilding matrix system: done." << endl;
        }
    }

    /**
     * @brief Generate summary at runtime.
     */
    void LinearElasticSystem::v_GenerateSummary(SolverUtils::SummaryList& s)
    {
        EquationSystem::SessionSummary(s);
    }

    /**
     * @brief Solve elliptic linear elastic system.
     *
     * The solve proceeds as follows:
     *
     * - Create a MultiRegions::GlobalLinSys object.
     * - Evaluate a forcing function from the session file.
     * - If the temperature term is enabled, evaluate this and add to forcing.
     * - Apply Dirichlet boundary conditions.
     * - Scatter forcing into the correct ordering according to the coupled
     *   assembly map.
     * - Do the solve.
     * - Scatter solution back to fields and backwards transform to physical
     *   space.
     */
    void LinearElasticSystem::v_DoSolve()
    {
        int i, j, k, l, nv;
        const int nVel = m_fields[0]->GetCoordim(0);

        // Build initial matrix system.
        BuildMatrixSystem();

        // Now we've got the matrix system set up, create a GlobalLinSys object.
        MultiRegions::GlobalLinSysKey key(
            StdRegions::eLinearAdvectionReaction, m_assemblyMap);        
        MultiRegions::GlobalLinSysSharedPtr linSys;

        // Currently either direct or iterative static condensation is
        // supported.
        if (m_assemblyMap->GetGlobalSysSolnType() == MultiRegions::eDirectStaticCond)
        {
            linSys = MemoryManager<
                MultiRegions::GlobalLinSysDirectStaticCond>::AllocateSharedPtr(
                    key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
                    m_assemblyMap);
        }
        else if (m_assemblyMap->GetGlobalSysSolnType() ==
                     MultiRegions::eIterativeStaticCond)
        {
            linSys = MemoryManager<
                MultiRegions::GlobalLinSysIterativeStaticCond>::AllocateSharedPtr(
                    key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
                    m_assemblyMap, MultiRegions::NullPreconditionerSharedPtr);
        }

        const int nCoeffs = m_fields[0]->GetNcoeffs();
        const int nGlobDofs = m_assemblyMap->GetNumGlobalCoeffs();

        //
        // -- Evaluate forcing functions
        //

        // Evaluate the forcing function from the XML file.
        Array<OneD, Array<OneD, NekDouble> > forcing(nVel);
        EvaluateFunction(forcing, "Forcing");

        // Add temperature term
        string tempEval;
        m_session->LoadSolverInfo("Temperature", tempEval, "None");
        if (tempEval == "Jacobian")
        {
            // Allocate storage
            m_temperature = Array<OneD, Array<OneD, NekDouble> >(nVel);

            for (nv = 0; nv < nVel; ++nv)
            {
                Array<OneD, NekDouble> tmp;
                m_temperature[nv] = Array<OneD, NekDouble>(
                    m_fields[nv]->GetNpoints());
                
                for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
                {
                    // Calculate element area
                    LocalRegions::ExpansionSharedPtr exp =
                        m_fields[0]->GetExp(i);
                    LibUtilities::PointsKeyVector pkey =
                        exp->GetPointsKeys();
                    Array<OneD, NekDouble> jac =
                        exp->GetMetricInfo()->GetJac(pkey);
                    
                    int offset = m_fields[0]->GetPhys_Offset(i);
                    Vmath::Smul(exp->GetTotPoints(), m_beta, jac, 1,
                                tmp = m_temperature[nv] + offset, 1);
                }
                m_fields[nv]->PhysDeriv(nv, m_temperature[nv], forcing[nv]);
            }
        }
        else if (tempEval == "Metric")
        {
            m_temperature = Array<OneD, Array<OneD, NekDouble> >(nVel);
            Array<TwoD, Array<OneD, NekDouble> > beta(nVel, nVel);

            for (nv = 0; nv < nVel; ++nv)
            {
                m_temperature[nv] = Array<OneD, NekDouble>(
                    m_fields[nv]->GetNpoints());
            }

            for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
            {
                // Calculate element area
                LocalRegions::ExpansionSharedPtr exp =
                    m_fields[0]->GetExp(i);
                LibUtilities::PointsKeyVector pkey = exp->GetPointsKeys();
                Array<OneD, NekDouble> jac =
                    exp->GetMetricInfo()->GetJac(pkey);
                Array<OneD, Array<OneD, Array<OneD, NekDouble> > > deriv
                    = exp->GetMetricInfo()->GetDeriv(pkey);
                int offset = m_fields[0]->GetPhys_Offset(i);
                
                Array<OneD, NekDouble> tmp_sv(jac.num_elements(),0.0);
                // Compute metric tensor
                Array<OneD, NekDouble> tmp(nVel*nVel, 0.0);
                for (j = 0; j < exp->GetTotPoints(); ++j)
                {
                    // TODO: fix this!
                    /*
                    for (k = 0; k < nVel; ++k)
                    {
                        for (l = 0; l < nVel; ++l)
                        {
                            tmp[k*nVel+l] = 0.0;
                            for (int m = 0; m < nVel; ++m)
                            {
                                tmp[k*nVel+l] += deriv[m][k][j] * deriv[l][m][j];
                            }
                            cout << tmp[k*nVel+l] << endl;
                        }
                    }*/
                    
                    tmp[0] = deriv[0][0][j] * deriv[0][0][j] +
                             deriv[0][1][j] * deriv[0][1][j];
                    tmp[1] = deriv[1][0][j] * deriv[0][0][j] +
                             deriv[1][1][j] * deriv[0][1][j];
                    tmp[2] = deriv[0][0][j] * deriv[1][0][j] +
                             deriv[0][1][j] * deriv[1][1][j];
                    tmp[3] = deriv[1][0][j] * deriv[1][0][j] +
                             deriv[1][1][j] * deriv[1][1][j];


                    // Compute eigenvalues/eigenvectors.
                    char jobvl = 'N', jobvr = 'V';
                    int worklen = 8*nVel, info;
                    Array<OneD, NekDouble> wi(nVel);
                    DNekMat evec(nVel, nVel, 0.0, eFULL);
                    DNekMat evecinv(nVel, nVel, 0.0, eFULL);
                    DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
                    Array<OneD, NekDouble> vl(nVel*nVel);
                    Array<OneD, NekDouble> work(worklen);
                    
                    //Lapack::Dgeev(jobvl, jobvr, nVel, &tmp[0], nVel,
                    //              &(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
                    //              &(evec.GetPtr())[0], nVel,
                    //              &work[0], worklen, info);

                    NekDouble T = -(tmp[0] + tmp[3]);
                    NekDouble D = tmp[0]*tmp[3] - tmp[1]*tmp[2];
                    
                    NekDouble Term = (T*T)-4*D;
                    eval(0,0) = (-T+sqrt(Term))/2;
                    eval(1,1) = (-T-sqrt(Term))/2;
                    
                    if (tmp[1] != 0)
                    {
                        evec(0,0) = tmp[1];
                        evec(1,0) = eval(0,0) - tmp[0];
                        
                        evec(0,1) = tmp[1];
                        evec(1,1) = eval(1,1) - tmp[0];
                    }
                    else if(tmp[2] != 0)
                    {
                        evec(0,0) = eval(0,0) - tmp[3];
                        evec(1,0) = tmp[2];
                        
                        evec(0,1) = eval(1,1) - tmp[3];
                        evec(1,1) = tmp[2];
                    }
                    else
                    {
                        evec(0,0) = 1.0;
                        evec(1,0) = 0.0;
                        
                        evec(0,1) = 0.0;
                        evec(1,1) = 1.0;
                    }
                    
                    NekDouble det = evec(0,0)*evec(1,1)-evec(0,1)*evec(1,0);
                    NekDouble fac = 1/det;
                    
                    // inverse of the eigenvectors
                    evecinv(0,0) =   fac*evec(1,1);
                    evecinv(0,1) =  -fac*evec(0,1);
                    evecinv(1,0) =  -fac*evec(1,0);
                    evecinv(1,1) =   fac*evec(0,0);
                    
                    eval(0,0) = m_beta * eval(0,0);
                    eval(1,1) = m_beta * eval(1,1);

                    // evecinv.Invert();
                    DNekMat beta = evec * eval * evecinv ;
                    
                    for (nv = 0; nv < nVel; ++nv)
                    {
                        m_temperature[nv][offset + j] =
                            (beta(nv,0) + beta(nv,1)) * jac[j];
                    }
                }
            }

            for (nv = 0; nv < nVel; ++nv)
            {
                m_fields[nv]->PhysDeriv(nv, m_temperature[nv], forcing[nv]);
            }
        }
        else if (tempEval != "None")
        {
            ASSERTL0(false, "Unknown temperature form: " + tempEval);
        }


        // Set up some temporary storage.
        //
        // - forCoeffs holds the forcing coefficients in a local ordering;
        //   however note that the ordering is different and dictated by the
        //   assembly map. We loop over each element, then the boundary degrees<commit_msg>Clean up metric temperature term<commit_after>///////////////////////////////////////////////////////////////////////////////
//
// File LinearElasticSystem.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: LinearElasticSystem solve routines 
//
///////////////////////////////////////////////////////////////////////////////

#include <LocalRegions/MatrixKey.h>
#include <MultiRegions/ContField2D.h>
#include <MultiRegions/ContField3D.h>
#include <MultiRegions/GlobalLinSysDirectStaticCond.h>
#include <MultiRegions/GlobalLinSysIterativeStaticCond.h>
#include <MultiRegions/Preconditioner.h>
#include <LinearElasticSolver/EquationSystems/LinearElasticSystem.h>

namespace Nektar
{
    string LinearElasticSystem::className = GetEquationSystemFactory().
        RegisterCreatorFunction("LinearElasticSystem",
                                LinearElasticSystem::create);

    LinearElasticSystem::LinearElasticSystem(
            const LibUtilities::SessionReaderSharedPtr& pSession)
        : EquationSystem(pSession)
    {
    }

    /**
     * @brief Set up the linear elasticity system.
     *
     * This routine loads the E and nu variables from the session file, creates
     * a coupled assembly map and creates containers for the statically
     * condensed block matrix system.
     */
    void LinearElasticSystem::v_InitObject()
    {
        EquationSystem::v_InitObject();

        const int nVel = m_fields[0]->GetCoordim(0);
        int n;

        ASSERTL0(nVel > 1, "Linear elastic solver not set up for"
                           " this dimension (only 2D/3D supported).");

        // Make sure that we have Young's modulus and Poisson ratio set.
        m_session->LoadParameter("E", m_E, 1.0);
        m_session->LoadParameter("nu", m_nu, 0.25);
        m_session->LoadParameter("Beta", m_beta, 1.0e-04);

        // Create a coupled assembly map which allows us to tie u, v and w
        // fields together.
        if (nVel == 2)
        {
            MultiRegions::ContField2DSharedPtr u = boost::dynamic_pointer_cast<
                MultiRegions::ContField2D>(m_fields[0]);
            m_assemblyMap = MemoryManager<CoupledAssemblyMap>
                ::AllocateSharedPtr(m_session,
                                    m_graph,
                                    u->GetLocalToGlobalMap(),
                                    m_boundaryConditions,
                                    m_fields);
        }
        
        if (nVel == 3)
        {
            MultiRegions::ContField3DSharedPtr u = boost::dynamic_pointer_cast<
                MultiRegions::ContField3D>(m_fields[0]);
            m_assemblyMap = MemoryManager<CoupledAssemblyMap>
                ::AllocateSharedPtr(m_session,
                                    m_graph,
                                    u->GetLocalToGlobalMap(),
                                    m_boundaryConditions,
                                    m_fields);
        }
        
        // Figure out size of our new matrix systems by looping over all
        // expansions and multiply number of coefficients by velocity
        // components.
        const int nEl = m_fields[0]->GetExpSize();
        LocalRegions::ExpansionSharedPtr exp;

        Array<OneD, unsigned int> sizeBnd(nEl);
        Array<OneD, unsigned int> sizeInt(nEl);

        // Allocate storage for boundary and interior map caches.
        m_bmap = Array<OneD, Array<OneD, unsigned int> >(nEl);
        m_imap = Array<OneD, Array<OneD, unsigned int> >(nEl);

        // Cache interior and boundary maps to avoid recalculating
        // constantly. Should really be handled by manager in StdRegions but not
        // really working yet.
        for (n = 0; n < nEl; ++n)
        {
            exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
            sizeBnd[n] = nVel * exp->NumBndryCoeffs();
            sizeInt[n] = nVel * exp->GetNcoeffs() - sizeBnd[n];
            exp->GetBoundaryMap(m_bmap[n]);
            exp->GetInteriorMap(m_imap[n]);
        }

        // Create block matrix storage for the statically condensed system.
        MatrixStorage blkmatStorage = eDIAGONAL;
        m_schurCompl = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            sizeBnd, sizeBnd, blkmatStorage);
        m_BinvD      = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            sizeBnd, sizeInt, blkmatStorage);
        m_C          = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            sizeInt, sizeBnd, blkmatStorage);
        m_Dinv       = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            sizeInt, sizeInt, blkmatStorage);
    }

    /**
     * @brief Build matrix system for linear elasticity equations.
     *
     * This routine constructs the matrix discretisation arising from the weak
     * formulation of the linear elasticity equations. The resulting matrix
     * system is then passed to LinearElasticSystem::SetStaticCondBlock in order
     * to construct the statically condensed system.
     *
     * All of the matrices involved in the construction of the divergence of the
     * stress tensor are Laplacian-like. We use the variable coefficient
     * functionality within the library to multiply the Laplacian by the
     * appropriate constants for all diagonal terms, and off-diagonal terms use
     * LinearElasticSystem::BuildLaplacianIJMatrix to construct matrices
     * containing
     *
     * \f[ \int \partial_{x_i} \phi_k \partial_{x_j} \phi_l \f]
     *
     * where mixed derivative terms are present. Symmetry (in terms of k,l) is
     * exploited to avoid constructing this matrix repeatedly.
     *
     * @todo Make static condensation optional and construct full system
     *       instead.
     */
    void LinearElasticSystem::BuildMatrixSystem()
    {
        const int nEl = m_fields[0]->GetExpSize();
        const int nVel = m_fields[0]->GetCoordim(0);
        
        LocalRegions::ExpansionSharedPtr exp;
        int n;

        // Factors map for matrix keys.
        StdRegions::ConstFactorMap factors;

        // Calculate various constants

        // mu
        NekDouble a = m_E * (1.0 - m_nu) / (1.0 + m_nu) / (1.0 - 2.0*m_nu);
        // lambda
        NekDouble b = m_E * 0.5 / (1.0 + m_nu);
        // lambda + mu?
        NekDouble c = m_E * m_nu / (1.0 + m_nu) / (1.0 - 2.0*m_nu);

        bool verbose = m_session->DefinesCmdLineArgument("verbose");
        bool root    = m_session->GetComm()->GetRank() == 0;

        // Loop over each element and construct matrices.
        if (nVel == 2)
        {
            for (n = 0; n < nEl; ++n)
            {
                exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
                const int nPhys = exp->GetTotPoints();
                Array<OneD, NekDouble> aArr(nPhys, a);
                Array<OneD, NekDouble> bArr(nPhys, b);
                
                StdRegions::VarCoeffMap varcoeffA, varcoeffD;
                varcoeffA[StdRegions::eVarCoeffD00] = aArr;
                varcoeffA[StdRegions::eVarCoeffD11] = bArr;
                varcoeffD[StdRegions::eVarCoeffD00] = bArr;
                varcoeffD[StdRegions::eVarCoeffD11] = aArr;

                LocalRegions::MatrixKey matkeyA(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffA);
                LocalRegions::MatrixKey matkeyD(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffD);

                /*
                 * mat holds the linear operator [ A B ] acting on [ u ].
                 *                               [ C D ]           [ v ]
                 */
                Array<TwoD, DNekMatSharedPtr> mat(2,2);
                mat[0][0] = exp->GenMatrix(matkeyA);
                mat[0][1] = BuildLaplacianIJMatrix(1, 0, c, exp);
                mat[1][0] = mat[0][1];
                mat[1][1] = exp->GenMatrix(matkeyD);

                // Set up the statically condensed block for this element.
                SetStaticCondBlock(n, exp, mat);

                if (verbose && root)
                {
                    cout << "\rBuilding matrix system: "
                         << (int)(100.0 * n / nEl) << "%" << flush;
                }
            }
        }
        else if (nVel == 3)
        {
            for (n = 0; n < nEl; ++n)
            {
                exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
                const int nPhys = exp->GetTotPoints();
                Array<OneD, NekDouble> aArr(nPhys, a);
                Array<OneD, NekDouble> bArr(nPhys, b);

                StdRegions::VarCoeffMap varcoeffA, varcoeffE, varcoeffI;
                varcoeffA[StdRegions::eVarCoeffD00] = aArr;
                varcoeffA[StdRegions::eVarCoeffD11] = bArr;
                varcoeffA[StdRegions::eVarCoeffD22] = bArr;
                varcoeffE[StdRegions::eVarCoeffD00] = bArr;
                varcoeffE[StdRegions::eVarCoeffD11] = aArr;
                varcoeffE[StdRegions::eVarCoeffD22] = bArr;
                varcoeffI[StdRegions::eVarCoeffD00] = bArr;
                varcoeffI[StdRegions::eVarCoeffD11] = bArr;
                varcoeffI[StdRegions::eVarCoeffD22] = aArr;

                LocalRegions::MatrixKey matkeyA(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffA);
                LocalRegions::MatrixKey matkeyE(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffE);
                LocalRegions::MatrixKey matkeyI(StdRegions::eLaplacian,
                                                exp->DetShapeType(),
                                                *exp, factors, varcoeffI);

                /*
                 * mat holds the linear operator [ A B C ] acting on [ u ].
                 *                               [ D E F ]           [ v ]
                 *                               [ G H I ]           [ w ]
                 */
                Array<TwoD, DNekMatSharedPtr> mat(3,3);
                mat[0][0] = exp->GenMatrix(matkeyA);
                mat[0][1] = BuildLaplacianIJMatrix(1, 0, c, exp);
                mat[0][2] = BuildLaplacianIJMatrix(2, 0, c, exp);
                
                mat[1][0] = mat[0][1];
                mat[1][1] = exp->GenMatrix(matkeyE);
                mat[1][2] = BuildLaplacianIJMatrix(2, 1, c, exp);
                
                mat[2][0] = mat[0][2];
                mat[2][1] = mat[1][2];
                mat[2][2] = exp->GenMatrix(matkeyI);
                
                // Set up the statically condensed block for this element.
                SetStaticCondBlock(n, exp, mat);

                if (verbose && root)
                {
                    cout << "\rBuilding matrix system: "
                         << (int)(100.0 * n / nEl) << "%" << flush;
                }
            }
        }

        if (verbose && root)
        {
            cout << "\rBuilding matrix system: done." << endl;
        }
    }

    /**
     * @brief Generate summary at runtime.
     */
    void LinearElasticSystem::v_GenerateSummary(SolverUtils::SummaryList& s)
    {
        EquationSystem::SessionSummary(s);
    }

    /**
     * @brief Solve elliptic linear elastic system.
     *
     * The solve proceeds as follows:
     *
     * - Create a MultiRegions::GlobalLinSys object.
     * - Evaluate a forcing function from the session file.
     * - If the temperature term is enabled, evaluate this and add to forcing.
     * - Apply Dirichlet boundary conditions.
     * - Scatter forcing into the correct ordering according to the coupled
     *   assembly map.
     * - Do the solve.
     * - Scatter solution back to fields and backwards transform to physical
     *   space.
     */
    void LinearElasticSystem::v_DoSolve()
    {
        int i, j, k, l, m, nv;
        const int nVel = m_fields[0]->GetCoordim(0);

        // Build initial matrix system.
        BuildMatrixSystem();

        // Now we've got the matrix system set up, create a GlobalLinSys object.
        MultiRegions::GlobalLinSysKey key(
            StdRegions::eLinearAdvectionReaction, m_assemblyMap);        
        MultiRegions::GlobalLinSysSharedPtr linSys;

        // Currently either direct or iterative static condensation is
        // supported.
        if (m_assemblyMap->GetGlobalSysSolnType() == MultiRegions::eDirectStaticCond)
        {
            linSys = MemoryManager<
                MultiRegions::GlobalLinSysDirectStaticCond>::AllocateSharedPtr(
                    key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
                    m_assemblyMap);
        }
        else if (m_assemblyMap->GetGlobalSysSolnType() ==
                     MultiRegions::eIterativeStaticCond)
        {
            linSys = MemoryManager<
                MultiRegions::GlobalLinSysIterativeStaticCond>::AllocateSharedPtr(
                    key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
                    m_assemblyMap, MultiRegions::NullPreconditionerSharedPtr);
        }

        const int nCoeffs = m_fields[0]->GetNcoeffs();
        const int nGlobDofs = m_assemblyMap->GetNumGlobalCoeffs();

        //
        // -- Evaluate forcing functions
        //

        // Evaluate the forcing function from the XML file.
        Array<OneD, Array<OneD, NekDouble> > forcing(nVel);
        EvaluateFunction(forcing, "Forcing");

        // Add temperature term
        string tempEval;
        m_session->LoadSolverInfo("Temperature", tempEval, "None");
        if (tempEval == "Jacobian")
        {
            // Allocate storage
            m_temperature = Array<OneD, Array<OneD, NekDouble> >(nVel);

            for (nv = 0; nv < nVel; ++nv)
            {
                Array<OneD, NekDouble> tmp;
                m_temperature[nv] = Array<OneD, NekDouble>(
                    m_fields[nv]->GetNpoints());
                
                for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
                {
                    // Calculate element area
                    LocalRegions::ExpansionSharedPtr exp =
                        m_fields[0]->GetExp(i);
                    LibUtilities::PointsKeyVector pkey =
                        exp->GetPointsKeys();
                    Array<OneD, NekDouble> jac =
                        exp->GetMetricInfo()->GetJac(pkey);
                    
                    int offset = m_fields[0]->GetPhys_Offset(i);
                    Vmath::Smul(exp->GetTotPoints(), m_beta, jac, 1,
                                tmp = m_temperature[nv] + offset, 1);
                }
                m_fields[nv]->PhysDeriv(nv, m_temperature[nv], forcing[nv]);
            }
        }
        else if (tempEval == "Metric")
        {
            m_temperature = Array<OneD, Array<OneD, NekDouble> >(nVel);
            Array<TwoD, Array<OneD, NekDouble> > beta(nVel, nVel);

            for (nv = 0; nv < nVel; ++nv)
            {
                m_temperature[nv] = Array<OneD, NekDouble>(
                    m_fields[nv]->GetNpoints());
            }

            for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
            {
                // Calculate element area
                LocalRegions::ExpansionSharedPtr exp =
                    m_fields[0]->GetExp(i);
                LibUtilities::PointsKeyVector pkey = exp->GetPointsKeys();
                Array<OneD, NekDouble> jac =
                    exp->GetMetricInfo()->GetJac(pkey);
                Array<OneD, Array<OneD, Array<OneD, NekDouble> > > deriv
                    = exp->GetMetricInfo()->GetDeriv(pkey);
                int offset = m_fields[0]->GetPhys_Offset(i);
                
                Array<OneD, NekDouble> tmp_sv(jac.num_elements(), 0.0);
                // Compute metric tensor
                Array<OneD, NekDouble> tmp(nVel*nVel, 0.0);
                for (j = 0; j < exp->GetTotPoints(); ++j)
                {
                    // Setting up the metric tensor
                    for (k = 0; k < nVel; ++k)
                    {
                        for (l = 0; l < nVel; ++l)
                        {
                            tmp[k*nVel+l] = 0.0;
                            
                            for (int m = 0; m < nVel; ++m)
                            {
                                tmp[k*nVel+l] += deriv[k][m][j] * deriv[l][m][j];
                            }
                        }
                    }
                
                    // Compute eigenvalues/eigenvectors.
                    char jobvl = 'N', jobvr = 'V';
                    int worklen = 8*nVel, info;
                    
                    DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
                    DNekMat evec(nVel, nVel, 0.0, eFULL);
                    DNekMat evecinv(nVel, nVel, 0.0, eFULL);
                    Array<OneD, NekDouble> vl(nVel*nVel);
                    Array<OneD, NekDouble> work(worklen);
                    Array<OneD, NekDouble> wi(nVel);
                    
                    Lapack::Dgeev(jobvl, jobvr, nVel, &tmp[0], nVel,
                                  &(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
                                  &(evec.GetPtr())[0], nVel,
                                  &work[0], worklen, info);
                    
                    evecinv = evec;
                    evecinv.Invert();

                    // rescaling of the eigenvalues
                    for (nv = 0; nv < nVel; ++nv)
                    {
                        eval(nv,nv) = m_beta * eval(nv,nv);
                    }
                
                    DNekMat beta = evec * eval * evecinv ;
                    NekDouble term = 0.0;
                    
                    for (nv = 0; nv < nVel; ++nv)
                    {
                        term = 0.0;
                        for (m = 0; m < nVel; ++m)
                        {
                            term += beta(nv,m);
                        }
                        m_temperature[nv][offset + j] =
                            (term) * jac[j];
                    }
                }
            }

            for (nv = 0; nv < nVel; ++nv)
            {
                m_fields[nv]->PhysDeriv(nv, m_temperature[nv], forcing[nv]);
            }
        }
        else if (tempEval != "None")
        {
            ASSERTL0(false, "Unknown temperature form: " + tempEval);
        }


        // Set up some temporary storage.
        //
        // - forCoeffs holds the forcing coefficients in a local ordering;
        //   however note that the ordering is different and dictated by the
        //   assembly map. We loop over each element, then the boundary degrees
{'fraction_non_alphanumeric': 0.11101994889706744, 'fraction_numerical': 0.01204904006188612, 'mean_word_length': 1.4489092996555684, 'pattern_counts': {'":': 0, '<': 175, '<?xml version=': 0, '>': 225, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 82, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '418583', 'n_tokens_mistral': 12450, 'n_tokens_neox': 11538, 'n_words': 2867}