<commit_before>// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

#include "../general/forall.hpp"
#include "bilininteg.hpp"
#include "gridfunc.hpp"

using namespace std;

namespace mfem
{
// PA DG Trace Integrator

// PA DG Trace Assemble 2D kernel for constant velocity
static void PADGTraceSetup2D(const int Q1D,
                             const int NF,
                             const Array<double> &w,
                             const Vector &det,
                             const Vector &nor,
                             const double rho,
                             const Vector &vel,
                             const double alpha,
                             const double beta,
                             Vector &op)
{
   const int VDIM = 2;

   auto d = Reshape(det.Read(), Q1D, NF);
   auto n = Reshape(nor.Read(), Q1D, VDIM, NF);
   // auto v = Reshape(vel.Read(), Q1D, VDIM, NF);//Not for constant velocity
   auto c = vel.Read();
   auto W = w.Read();
   auto qd = Reshape(op.Write(), Q1D, 2, 2, NF);

   MFEM_FORALL(f, NF,//can be optimized with Q1D thread for NF blocks
   {
      for (int q = 0; q < Q1D; ++q)
      {
         // const double dot = n(q,0,f) * c(q,0,f) + n(q,1,f) * c(q,1,f);
         const double dot = n(q,0,f) * c[0] + n(q,1,f) * c[1];
         // const double dot = c[0] + c[1];//just for debugging
         const double abs = dot > 0.0 ? dot : -dot;
         const double w = W[q]*d(q,f);
         qd(q,0,0,f) = w*( alpha/2 * dot + beta * abs );
         qd(q,1,0,f) = w*( alpha/2 * dot - beta * abs );
         qd(q,0,1,f) = w*(-alpha/2 * dot - beta * abs );
         qd(q,1,1,f) = w*(-alpha/2 * dot + beta * abs );
      }
   });
}

// PA DG Trace Assemble 2D kernel for constant velocity
static void PADGTraceSetup3D(const int Q1D,
                             const int NF,
                             const Array<double> &w,
                             const Vector &det,
                             const Vector &nor,
                             const double rho,
                             const Vector &vel,
                             const double alpha,
                             const double beta,
                             Vector &op)
{
   const int VDIM = 3;

   auto d = Reshape(det.Read(), Q1D, Q1D, NF);
   auto n = Reshape(nor.Read(), Q1D, Q1D, VDIM, NF);
   // auto c = Reshape(vel.Read(), Q1D, Q1D, VDIM, NF);
   auto c = vel.Read();
   auto W = w.Read();
   auto qd = Reshape(op.Write(), Q1D, Q1D, 2, 2, NF);

   MFEM_FORALL(f, NF,//can be optimized with Q1D*Q1D threads for NF blocks
   {
      for (int q1 = 0; q1 < Q1D; ++q1)
      {
         for (int q2 = 0; q2 < Q1D; ++q2)
         {
            // const double dot = n(q,0,f) * c(q,0,f) + n(q,1,f) * c(q,1,f) + n(q,2,f) * c(q,2,f);
            const double dot = n(q1,q2,0,f) * c[0] + n(q1,q2,1,f) * c[1] + n(q1,q2,2,f) * c[2];
            const double abs = dot > 0.0 ? dot : -dot;
            const double w = W[q1]*W[q2]*d(q1,q2,f);
            qd(q1,q2,0,0,f) = w*( alpha/2 * dot + beta * abs );
            qd(q1,q2,1,0,f) = w*( alpha/2 * dot - beta * abs );
            qd(q1,q2,0,1,f) = w*(-alpha/2 * dot - beta * abs );
            qd(q1,q2,1,1,f) = w*(-alpha/2 * dot + beta * abs );
         }
      }
   });
}

static void PADGTraceSetup(const int dim,
                           const int D1D,
                           const int Q1D,
                           const int NF,
                           const Array<double> &W,
                           const Vector &det,
                           const Vector &nor,
                           const double rho,
                           const Vector &u,
                           const double alpha,
                           const double beta,
                           Vector &op)
{
   if (dim == 1) { MFEM_ABORT("dim==1 not supported in PADGTraceSetup"); }
   if (dim == 2)
   {
      PADGTraceSetup2D(Q1D, NF, W, det, nor, rho, u, alpha, beta, op);
   }
   if (dim == 3)
   {
      PADGTraceSetup3D(Q1D, NF, W, det, nor, rho, u, alpha, beta, op);
   }
}

void DGTraceIntegrator::AssemblePA(const FiniteElementSpace &fes)
{
   // Assumes tensor-product elements
   Mesh *mesh = fes.GetMesh();
   const FiniteElement &el = *fes.GetFaceElement(0);
   const IntegrationRule *ir = &GetRule(el, el);//IntRule ? IntRule : &GetRule(el, el);
   const int dims = el.GetDim();
   const int symmDims = 4;
   const int nq = ir->GetNPoints();
   dim = mesh->Dimension();
   nf = mesh->GetNumFaces();
   geom = mesh->GetFaceGeometricFactors(*ir,
      FaceGeometricFactors::DETERMINANTS | FaceGeometricFactors::NORMALS);
   maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
   dofs1D = maps->ndof;
   quad1D = maps->nqpt;
   pa_data.SetSize(symmDims * nq * nf, Device::GetMemoryType());
   double r;<commit_msg>Workaround for GetFaceElement<commit_after>// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

#include "../general/forall.hpp"
#include "bilininteg.hpp"
#include "gridfunc.hpp"

using namespace std;

namespace mfem
{
// PA DG Trace Integrator

// PA DG Trace Assemble 2D kernel for constant velocity
static void PADGTraceSetup2D(const int Q1D,
                             const int NF,
                             const Array<double> &w,
                             const Vector &det,
                             const Vector &nor,
                             const double rho,
                             const Vector &vel,
                             const double alpha,
                             const double beta,
                             Vector &op)
{
   const int VDIM = 2;

   auto d = Reshape(det.Read(), Q1D, NF);
   auto n = Reshape(nor.Read(), Q1D, VDIM, NF);
   // auto v = Reshape(vel.Read(), Q1D, VDIM, NF);//Not for constant velocity
   auto c = vel.Read();
   auto W = w.Read();
   auto qd = Reshape(op.Write(), Q1D, 2, 2, NF);

   MFEM_FORALL(f, NF,//can be optimized with Q1D thread for NF blocks
   {
      for (int q = 0; q < Q1D; ++q)
      {
         // const double dot = n(q,0,f) * c(q,0,f) + n(q,1,f) * c(q,1,f);
         const double dot = n(q,0,f) * c[0] + n(q,1,f) * c[1];
         // const double dot = c[0] + c[1];//just for debugging
         const double abs = dot > 0.0 ? dot : -dot;
         const double w = W[q]*d(q,f);
         qd(q,0,0,f) = w*( alpha/2 * dot + beta * abs );
         qd(q,1,0,f) = w*( alpha/2 * dot - beta * abs );
         qd(q,0,1,f) = w*(-alpha/2 * dot - beta * abs );
         qd(q,1,1,f) = w*(-alpha/2 * dot + beta * abs );
      }
   });
}

// PA DG Trace Assemble 2D kernel for constant velocity
static void PADGTraceSetup3D(const int Q1D,
                             const int NF,
                             const Array<double> &w,
                             const Vector &det,
                             const Vector &nor,
                             const double rho,
                             const Vector &vel,
                             const double alpha,
                             const double beta,
                             Vector &op)
{
   const int VDIM = 3;

   auto d = Reshape(det.Read(), Q1D, Q1D, NF);
   auto n = Reshape(nor.Read(), Q1D, Q1D, VDIM, NF);
   // auto c = Reshape(vel.Read(), Q1D, Q1D, VDIM, NF);
   auto c = vel.Read();
   auto W = w.Read();
   auto qd = Reshape(op.Write(), Q1D, Q1D, 2, 2, NF);

   MFEM_FORALL(f, NF,//can be optimized with Q1D*Q1D threads for NF blocks
   {
      for (int q1 = 0; q1 < Q1D; ++q1)
      {
         for (int q2 = 0; q2 < Q1D; ++q2)
         {
            // const double dot = n(q,0,f) * c(q,0,f) + n(q,1,f) * c(q,1,f) + n(q,2,f) * c(q,2,f);
            const double dot = n(q1,q2,0,f) * c[0] + n(q1,q2,1,f) * c[1] + n(q1,q2,2,f) * c[2];
            const double abs = dot > 0.0 ? dot : -dot;
            const double w = W[q1]*W[q2]*d(q1,q2,f);
            qd(q1,q2,0,0,f) = w*( alpha/2 * dot + beta * abs );
            qd(q1,q2,1,0,f) = w*( alpha/2 * dot - beta * abs );
            qd(q1,q2,0,1,f) = w*(-alpha/2 * dot - beta * abs );
            qd(q1,q2,1,1,f) = w*(-alpha/2 * dot + beta * abs );
         }
      }
   });
}

static void PADGTraceSetup(const int dim,
                           const int D1D,
                           const int Q1D,
                           const int NF,
                           const Array<double> &W,
                           const Vector &det,
                           const Vector &nor,
                           const double rho,
                           const Vector &u,
                           const double alpha,
                           const double beta,
                           Vector &op)
{
   if (dim == 1) { MFEM_ABORT("dim==1 not supported in PADGTraceSetup"); }
   if (dim == 2)
   {
      PADGTraceSetup2D(Q1D, NF, W, det, nor, rho, u, alpha, beta, op);
   }
   if (dim == 3)
   {
      PADGTraceSetup3D(Q1D, NF, W, det, nor, rho, u, alpha, beta, op);
   }
}

void DGTraceIntegrator::AssemblePA(const FiniteElementSpace &fes)
{
   // Assumes tensor-product elements
   Mesh *mesh = fes.GetMesh();
   const FiniteElement &el = *fes.GetFaceElement(0);
   const IntegrationRule *ir = &GetRule(el, el);//IntRule ? IntRule : &GetRule(el, el);
   // Note: GetFE(0) and GetFaceElement(0) return elements with *different*
   // basis types (e.g. GaussLobatto and GaussLegendre)
   // TODO: find a better solution for this
   const FiniteElement &el_tmp = *fes.GetFE(0);
   const IntegrationRule &ir_tmp = IntRules.Get(el_tmp.GetGeomType(), ir->GetOrder());
   const int dims = el.GetDim();
   const int symmDims = 4;
   const int nq = ir->GetNPoints();
   dim = mesh->Dimension();
   nf = mesh->GetNumFaces();
   geom = mesh->GetFaceGeometricFactors(*ir,
      FaceGeometricFactors::DETERMINANTS | FaceGeometricFactors::NORMALS);
   maps = &el_tmp.GetDofToQuad(ir_tmp, DofToQuad::TENSOR);
   dofs1D = maps->ndof;
   quad1D = maps->nqpt;
   pa_data.SetSize(symmDims * nq * nf, Device::GetMemoryType());
   double r;
{'fraction_non_alphanumeric': 0.13464004396409598, 'fraction_numerical': 0.025370947059901082, 'mean_word_length': 1.6887466141344496, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1750390', 'n_tokens_mistral': 3879, 'n_tokens_neox': 3731, 'n_words': 1065}