<commit_before>//                               ETHOS Example 1
//
// Compile with: make exsm
//
// Sample runs:  exsm -m tipton.mesh
//
// Description: This example code performs a simple mesh smoothing based on a
//              topologically defined "mesh Laplacian" matrix.
//
//              The example highlights meshes with curved elements, the
//              assembling of a custom finite element matrix, the use of vector
//              finite element spaces, the definition of different spaces and
//              grid functions on the same mesh, and the setting of values by
//              iterating over the interior and the boundary elements.

#include "mfem.hpp"
#include <fstream>
#include <iostream>
#include <limits>
#include <unistd.h>

using namespace mfem;

using namespace std;

// This is where K10 methods start
// Implement Relaxed Newton solver to avoid mesh becoming invalid during the
// optimization process
class RelaxedNewtonSolver : public IterativeSolver
{
protected:
    mutable Vector r, c;
public:
    RelaxedNewtonSolver() { }
    
#ifdef MFEM_USE_MPI
    RelaxedNewtonSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
#endif
        virtual void SetOperator(const Operator &op);
    
        virtual void SetSolver(Solver &solver) { prec = &solver; }
    
        virtual void Mult(const Vector &b, Vector &x) const;
    
        virtual void Mult2(const Vector &b, Vector &x, ParMesh &pmesh,
                           const IntegrationRule &ir, int *itnums,
                           const ParNonlinearForm &nlf) const;
};

void RelaxedNewtonSolver::SetOperator(const Operator &op)
{
    oper = &op;
    height = op.Height();
    width = op.Width();
    MFEM_ASSERT(height == width, "square Operator is required.");
    
    r.SetSize(width);
    c.SetSize(width);
}

void RelaxedNewtonSolver::Mult(const Vector &b, Vector &x) const
{
    return;
}

void RelaxedNewtonSolver::Mult2(const Vector &b, Vector &x,
                                ParMesh &pmesh, const IntegrationRule &ir,
                                int *itnums, const ParNonlinearForm &nlf) const
{
    int myid;
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    
    MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
    MFEM_ASSERT(prec != NULL, "the Solver is not set (use SetSolver).");
    
    int it;
    double norm, norm_goal;
    bool have_b = (b.Size() == Height());
    
    if (!iterative_mode)
    {
        x = 0.0;
    }
    
    oper->Mult(x, r);
    if (have_b)
    {
        r -= b;
    }
    
    norm = Norm(r);
    norm_goal = std::max(rel_tol*norm, abs_tol);
    
    prec->iterative_mode = false;
    
    // x_{i+1} = x_i - [DF(x_i)]^{-1} [F(x_i)-b]
    for (it = 0; true; it++)
    {
        MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
        if (print_level >= 0)
        {
            cout << "Newton iteration " << setw(2) << it
                 << " : ||r|| = " << norm << '\n';
        }
        
        if (norm <= norm_goal)
        {
            converged = 1;
            break;
        }
        
        if (it >= max_iter)
        {
            converged = 0;
            break;
        }
        
        prec->SetOperator(oper->GetGradient(x));
        prec->Mult(r, c);  // c = [DF(x_i)]^{-1} [F(x_i)-b]
        
        //k10 some changes
        const int NE = pmesh.GetNE();
    
        Vector xsav = x; //create a copy of x
        Vector csav = c;
        int tchk = 0;
        int iters = 0;
        double alpha = 1.;
        int jachk = 1;
        double initenergy = 0.0;
        double finenergy = 0.0;
        int nanchk;
        
        initenergy = nlf.GetEnergy(x);
        
        while (tchk !=1 && iters < 20)
        {
            iters += 1;
            //cout << "number of iters is " << iters << "\n";
            jachk = 1;
            c = csav;
            x = xsav;
            add (xsav,-alpha,csav,x);

            ParGridFunction x_gf(nlf.ParFESpace());
            x_gf.Distribute(x);

            finenergy = nlf.GetEnergy(x);
            nanchk = isnan(finenergy);
            for (int i = 0; i < NE; i++)
            {
                const FiniteElement &fe = *x_gf.ParFESpace()->GetFE(i);
                const int dim = fe.GetDim(), nsp = ir.GetNPoints(),
                dof = fe.GetDof();
                
                DenseMatrix Jtr(dim);
                
                DenseMatrix dshape(dof, dim), pos(dof, dim);
                Array<int> xdofs(dof * dim);
                Vector posV(pos.Data(), dof * dim);
                
                x_gf.ParFESpace()->GetElementVDofs(i, xdofs);
                x_gf.GetSubVector(xdofs, posV);
                for (int j = 0; j < nsp; j++)
                {
                    fe.CalcDShape(ir.IntPoint(j), dshape);
                    MultAtB(pos, dshape, Jtr);
                    double det = Jtr.Det();
                }
            }
            int jachkglob;
            MPI_Allreduce(&jachk, &jachkglob, 1, MPI_INT, MPI_MIN,
                          pmesh.GetComm());
            int nanchkglob;
            MPI_Allreduce(&nanchk, &nanchkglob, 1, MPI_INT, MPI_MIN,
                          pmesh.GetComm());
            
            //finenergy = initenergy - 1; //WARNING: THIS IS JUST FOR TIPTONUNI.MESH
            if (finenergy>1.0*initenergy || nanchkglob!=0 || jachkglob==0)
            {
                tchk = 0;
                alpha *= 0.5;
            }
            else
            {
                tchk = 1;<commit_msg>Fixed a bug.<commit_after>//                               ETHOS Example 1
//
// Compile with: make exsm
//
// Sample runs:  exsm -m tipton.mesh
//
// Description: This example code performs a simple mesh smoothing based on a
//              topologically defined "mesh Laplacian" matrix.
//
//              The example highlights meshes with curved elements, the
//              assembling of a custom finite element matrix, the use of vector
//              finite element spaces, the definition of different spaces and
//              grid functions on the same mesh, and the setting of values by
//              iterating over the interior and the boundary elements.

#include "mfem.hpp"
#include <fstream>
#include <iostream>
#include <limits>
#include <unistd.h>

using namespace mfem;

using namespace std;

// This is where K10 methods start
// Implement Relaxed Newton solver to avoid mesh becoming invalid during the
// optimization process
class RelaxedNewtonSolver : public IterativeSolver
{
protected:
    mutable Vector r, c;
public:
    RelaxedNewtonSolver() { }
    
#ifdef MFEM_USE_MPI
    RelaxedNewtonSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
#endif
        virtual void SetOperator(const Operator &op);
    
        virtual void SetSolver(Solver &solver) { prec = &solver; }
    
        virtual void Mult(const Vector &b, Vector &x) const;
    
        virtual void Mult2(const Vector &b, Vector &x, ParMesh &pmesh,
                           const IntegrationRule &ir, int *itnums,
                           const ParNonlinearForm &nlf) const;
};

void RelaxedNewtonSolver::SetOperator(const Operator &op)
{
    oper = &op;
    height = op.Height();
    width = op.Width();
    MFEM_ASSERT(height == width, "square Operator is required.");
    
    r.SetSize(width);
    c.SetSize(width);
}

void RelaxedNewtonSolver::Mult(const Vector &b, Vector &x) const
{
    return;
}

void RelaxedNewtonSolver::Mult2(const Vector &b, Vector &x,
                                ParMesh &pmesh, const IntegrationRule &ir,
                                int *itnums, const ParNonlinearForm &nlf) const
{
    int myid;
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    
    MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
    MFEM_ASSERT(prec != NULL, "the Solver is not set (use SetSolver).");
    
    int it;
    double norm, norm_goal;
    bool have_b = (b.Size() == Height());
    
    if (!iterative_mode)
    {
        x = 0.0;
    }
    
    oper->Mult(x, r);
    if (have_b)
    {
        r -= b;
    }
    
    norm = Norm(r);
    norm_goal = std::max(rel_tol*norm, abs_tol);
    
    prec->iterative_mode = false;
    
    // x_{i+1} = x_i - [DF(x_i)]^{-1} [F(x_i)-b]
    for (it = 0; true; it++)
    {
        MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
        if (print_level >= 0)
        {
            cout << "Newton iteration " << setw(2) << it
                 << " : ||r|| = " << norm << '\n';
        }
        
        if (norm <= norm_goal)
        {
            converged = 1;
            break;
        }
        
        if (it >= max_iter)
        {
            converged = 0;
            break;
        }
        
        prec->SetOperator(oper->GetGradient(x));
        prec->Mult(r, c);  // c = [DF(x_i)]^{-1} [F(x_i)-b]
        
        //k10 some changes
        const int NE = pmesh.GetNE();
    
        Vector xsav = x; //create a copy of x
        Vector csav = c;
        int tchk = 0;
        int iters = 0;
        double alpha = 1.;
        int jachk = 1;
        double initenergy = 0.0;
        double finenergy = 0.0;
        int nanchk = 0;
        
        initenergy = nlf.GetEnergy(x);
        
        while (tchk !=1 && iters < 20)
        {
            iters += 1;
            //cout << "number of iters is " << iters << "\n";
            jachk = 1;
            c = csav;
            x = xsav;
            add (xsav,-alpha,csav,x);

            ParGridFunction x_gf(nlf.ParFESpace());
            x_gf.Distribute(x);

            finenergy = nlf.GetEnergy(x);
            nanchk = isnan(finenergy);
            for (int i = 0; i < NE; i++)
            {
                const FiniteElement &fe = *x_gf.ParFESpace()->GetFE(i);
                const int dim = fe.GetDim(), nsp = ir.GetNPoints(),
                dof = fe.GetDof();
                
                DenseMatrix Jtr(dim);
                
                DenseMatrix dshape(dof, dim), pos(dof, dim);
                Array<int> xdofs(dof * dim);
                Vector posV(pos.Data(), dof * dim);
                
                x_gf.ParFESpace()->GetElementVDofs(i, xdofs);
                x_gf.GetSubVector(xdofs, posV);
                for (int j = 0; j < nsp; j++)
                {
                    fe.CalcDShape(ir.IntPoint(j), dshape);
                    MultAtB(pos, dshape, Jtr);
                    double det = Jtr.Det();
                    if (det<=0.)
                    {
                        jachk = 0;
                    }
                }
            }
            int jachkglob;
            MPI_Allreduce(&jachk, &jachkglob, 1, MPI_INT, MPI_MIN,
                          pmesh.GetComm());
            int nanchkglob;
            MPI_Allreduce(&nanchk, &nanchkglob, 1, MPI_INT, MPI_MAX,
                          pmesh.GetComm());
            
            //finenergy = initenergy - 1; //WARNING: THIS IS JUST FOR TIPTONUNI.MESH
            if (finenergy>1.0*initenergy || nanchkglob!=0 || jachkglob==0)
            {
                tchk = 0;
                alpha *= 0.5;
            }
            else
            {
                tchk = 1;
{'fraction_non_alphanumeric': 0.11601615074024225, 'fraction_numerical': 0.007985643786451323, 'mean_word_length': 1.491282968261064, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5586742', 'n_tokens_mistral': 3623, 'n_tokens_neox': 3445, 'n_words': 918}