<commit_before>// Copyright (c) 2012-2016, <NAME>
// All rights reserved.
// 
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "integrators/gauss_rules.h"
#include "polywog/mlpg_quadrature.h"

typedef struct
{
  point_cloud_t* cloud;
  real_t* extents;
  int N;
  real_t ratio; 
} mlpg_simple_t;

static mlpg_simple_t* mlpg_simple_new(point_cloud_t* cloud,
                                      real_t* extents,
                                      int num_points,
                                      real_t ratio)
{
  ASSERT(num_points > 0);
  ASSERT(ratio > 0.0);

  mlpg_simple_t* mlpg = polymec_malloc(sizeof(mlpg_simple_t));
  mlpg->cloud = cloud;
  mlpg->extents = extents;
  mlpg->N = num_points;
  mlpg->ratio = ratio;
  return mlpg;
}

static void mlpg_simple_free(void* context)
{
  polymec_free(context);
}

static int cube_surf_num_quad_points(void* context, int i)
{
  mlpg_simple_t* mlpg = context;
  return 6 * mlpg->N * mlpg->N;
}

static void cube_surf_get_quad(void* context, int i, point_t* points, real_t* weights, vector_t* normals)
{
  mlpg_simple_t* mlpg = context;
  real_t gauss_pts[mlpg->N], gauss_wts[mlpg->N];
  get_gauss_legendre_points(mlpg->N, gauss_wts, gauss_wts);
  point_t* xi = &mlpg->cloud->points[i];
  real_t hi = mlpg->extents[i];
  real_t L = mlpg->ratio * hi;
  real_t a = xi->x - 0.5*L, b = xi->x + 0.5*L;
  real_t c = xi->y - 0.5*L, d = xi->y + 0.5*L;
  real_t e = xi->z - 0.5*L, f = xi->z + 0.5*L;
  int N = mlpg->N;
  int m = 0;
  for (int jj = 0; jj < mlpg->N; ++jj)
  {
    for (int kk = 0; kk < mlpg->N; ++kk, ++m)
    {
      real_t xj = (b - a) * gauss_pts[jj] + a;
      real_t xk = (b - a) * gauss_pts[kk] + a;
      real_t yj = (d - c) * gauss_pts[jj] + c;
      real_t yk = (d - c) * gauss_pts[kk] + c;
      real_t zj = (f - e) * gauss_pts[jj] + e;
      real_t zk = (f - e) * gauss_pts[kk] + e;
      real_t wjk = gauss_wts[jj] * gauss_wts[kk];

      // X coordinates for each of the faces.
      points[m].x       = a;
      normals[m].x      = -1.0;
      points[m+N*N].x   = b;
      normals[m+N*N].x  = 1.0;
      points[m+2*N*N].x = xj;
      normals[m+2*N*N].x = 0.0;
      points[m+3*N*N].x = xj;
      normals[m+3*N*N].x = 0.0;
      points[m+4*N*N].x = xk;
      normals[m+4*N*N].x = 0.0;
      points[m+5*N*N].x = xk;
      normals[m+5*N*N].x = 0.0;

      // Y coordinates for each of the faces.
      points[m].y = yj;
      normals[m].y = 0.0;
      points[m+N*N].y = yj;
      normals[m+N*N].y = 0.0;
      points[m+2*N*N].y = c;
      normals[m+2*N*N].y = -1.0;
      points[m+3*N*N].y = d;
      normals[m+3*N*N].y = 1.0;
      points[m+4*N*N].y = yk;
      normals[m+4*N*N].y = 0.0;
      points[m+5*N*N].y = yk;
      normals[m+5*N*N].y = 0.0;

      // Z coordinates for each of the faces.
      points[m].z = zj;
      normals[m].z = 0.0;
      points[m+N*N].z = zj;
      normals[m+N*N].z = 0.0;
      points[m+2*N*N].z = zk;
      normals[m+2*N*N].z = 0.0;
      points[m+3*N*N].z = zk;
      normals[m+3*N*N].z = 0.0;
      points[m+4*N*N].z = e;
      normals[m+4*N*N].z = -1.0;
      points[m+5*N*N].z = f;
      normals[m+5*N*N].z = 1.0;

      // Weights.
      weights[m] = wjk;
      weights[m+N*N] = wjk;
      weights[m+2*N*N] = wjk;
      weights[m+3*N*N] = wjk;
      weights[m+4*N*N] = wjk;
      weights[m+5*N*N] = wjk;
    }
  }
}

surface_integral_t* mlpg_cube_surface_integral_new(point_cloud_t* cloud,
                                                   real_t* extents,
                                                   int num_points,
                                                   real_t side_to_extent_ratio)
{
  mlpg_simple_t* mlpg = mlpg_simple_new(cloud, extents, num_points, 
                                        side_to_extent_ratio);
  surface_integral_vtable vtable = {.num_quad_points = cube_surf_num_quad_points,
                                    .get_quadrature = cube_surf_get_quad,
                                    .dtor = mlpg_simple_free};
  char name[1025];
  snprintf(name, 1024, "MLPG cube surface integral (N = %d, side/extent = %g)", 
           num_points, side_to_extent_ratio);
  return surface_integral_new(name, mlpg, vtable);
}

static int cube_vol_num_quad_points(void* context, int i)
{
  mlpg_simple_t* mlpg = context;
  return mlpg->N * mlpg->N * mlpg->N;
}

static void cube_vol_get_quad(void* context, int i, point_t* points, real_t* weights)
{
  mlpg_simple_t* mlpg = context;
  real_t gauss_pts[mlpg->N], gauss_wts[mlpg->N];
  get_gauss_legendre_points(mlpg->N, gauss_wts, gauss_wts);
  point_t* xi = &mlpg->cloud->points[i];
  real_t hi = mlpg->extents[i];
  real_t L = mlpg->ratio * hi;
  real_t a = xi->x - 0.5*L, b = xi->x + 0.5*L;
  real_t c = xi->y - 0.5*L, d = xi->y + 0.5*L;
  real_t e = xi->z - 0.5*L, f = xi->z + 0.5*L;
  real_t V = (b - a) * (d - c) * (f - e);
  int m = 0;
  for (int ii = 0; ii < mlpg->N; ++ii)
  {
    for (int jj = 0; jj < mlpg->N; ++jj)
    {
      for (int kk = 0; kk < mlpg->N; ++kk, ++m)
      {
        points[m].x = (b - a) * gauss_pts[ii] + a;
        points[m].y = (d - c) * gauss_pts[jj] + c;
        points[m].z = (f - e) * gauss_pts[kk] + e;
        weights[m] = gauss_wts[ii]*gauss_wts[jj]*gauss_wts[kk] * V;
      }
    }
  }
}

volume_integral_t* mlpg_cube_volume_integral_new(point_cloud_t* cloud,
                                                 real_t* extents,
                                                 int num_points,
                                                 real_t side_to_extent_ratio)
{
  mlpg_simple_t* mlpg = mlpg_simple_new(cloud, extents, num_points, 
                                        side_to_extent_ratio);
  volume_integral_vtable vtable = {.num_quad_points = cube_vol_num_quad_points,<commit_msg>Fixed a bug in the MLPG quadrature class.<commit_after>// Copyright (c) 2012-2016, <NAME>
// All rights reserved.
// 
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "integrators/gauss_rules.h"
#include "polywog/mlpg_quadrature.h"

typedef struct
{
  point_cloud_t* cloud;
  real_t* extents;
  int N;
  real_t ratio; 
} mlpg_simple_t;

static mlpg_simple_t* mlpg_simple_new(point_cloud_t* cloud,
                                      real_t* extents,
                                      int num_points,
                                      real_t ratio)
{
  ASSERT(num_points > 0);
  ASSERT(ratio > 0.0);

  mlpg_simple_t* mlpg = polymec_malloc(sizeof(mlpg_simple_t));
  mlpg->cloud = cloud;
  mlpg->extents = extents;
  mlpg->N = num_points;
  mlpg->ratio = ratio;
  return mlpg;
}

static void mlpg_simple_free(void* context)
{
  polymec_free(context);
}

static int cube_surf_num_quad_points(void* context, int i)
{
  mlpg_simple_t* mlpg = context;
  return 6 * mlpg->N * mlpg->N;
}

static void cube_surf_get_quad(void* context, int i, point_t* points, real_t* weights, vector_t* normals)
{
  mlpg_simple_t* mlpg = context;
  real_t gauss_pts[mlpg->N], gauss_wts[mlpg->N];
  get_gauss_legendre_points(mlpg->N, gauss_pts, gauss_wts);
  point_t* xi = &mlpg->cloud->points[i];
  real_t hi = mlpg->extents[i];
  real_t L = mlpg->ratio * hi;
  real_t a = xi->x - 0.5*L, b = xi->x + 0.5*L;
  real_t c = xi->y - 0.5*L, d = xi->y + 0.5*L;
  real_t e = xi->z - 0.5*L, f = xi->z + 0.5*L;
  int N = mlpg->N;
  int m = 0;
  for (int jj = 0; jj < mlpg->N; ++jj)
  {
    for (int kk = 0; kk < mlpg->N; ++kk, ++m)
    {
      real_t xj = (b - a) * gauss_pts[jj] + a;
      real_t xk = (b - a) * gauss_pts[kk] + a;
      real_t yj = (d - c) * gauss_pts[jj] + c;
      real_t yk = (d - c) * gauss_pts[kk] + c;
      real_t zj = (f - e) * gauss_pts[jj] + e;
      real_t zk = (f - e) * gauss_pts[kk] + e;
      real_t wjk = gauss_wts[jj] * gauss_wts[kk];

      // X coordinates for each of the faces.
      points[m].x       = a;
      normals[m].x      = -1.0;
      points[m+N*N].x   = b;
      normals[m+N*N].x  = 1.0;
      points[m+2*N*N].x = xj;
      normals[m+2*N*N].x = 0.0;
      points[m+3*N*N].x = xj;
      normals[m+3*N*N].x = 0.0;
      points[m+4*N*N].x = xk;
      normals[m+4*N*N].x = 0.0;
      points[m+5*N*N].x = xk;
      normals[m+5*N*N].x = 0.0;

      // Y coordinates for each of the faces.
      points[m].y = yj;
      normals[m].y = 0.0;
      points[m+N*N].y = yj;
      normals[m+N*N].y = 0.0;
      points[m+2*N*N].y = c;
      normals[m+2*N*N].y = -1.0;
      points[m+3*N*N].y = d;
      normals[m+3*N*N].y = 1.0;
      points[m+4*N*N].y = yk;
      normals[m+4*N*N].y = 0.0;
      points[m+5*N*N].y = yk;
      normals[m+5*N*N].y = 0.0;

      // Z coordinates for each of the faces.
      points[m].z = zj;
      normals[m].z = 0.0;
      points[m+N*N].z = zj;
      normals[m+N*N].z = 0.0;
      points[m+2*N*N].z = zk;
      normals[m+2*N*N].z = 0.0;
      points[m+3*N*N].z = zk;
      normals[m+3*N*N].z = 0.0;
      points[m+4*N*N].z = e;
      normals[m+4*N*N].z = -1.0;
      points[m+5*N*N].z = f;
      normals[m+5*N*N].z = 1.0;

      // Weights.
      weights[m] = wjk;
      weights[m+N*N] = wjk;
      weights[m+2*N*N] = wjk;
      weights[m+3*N*N] = wjk;
      weights[m+4*N*N] = wjk;
      weights[m+5*N*N] = wjk;
    }
  }
}

surface_integral_t* mlpg_cube_surface_integral_new(point_cloud_t* cloud,
                                                   real_t* extents,
                                                   int num_points,
                                                   real_t side_to_extent_ratio)
{
  mlpg_simple_t* mlpg = mlpg_simple_new(cloud, extents, num_points, 
                                        side_to_extent_ratio);
  surface_integral_vtable vtable = {.num_quad_points = cube_surf_num_quad_points,
                                    .get_quadrature = cube_surf_get_quad,
                                    .dtor = mlpg_simple_free};
  char name[1025];
  snprintf(name, 1024, "MLPG cube surface integral (N = %d, side/extent = %g)", 
           num_points, side_to_extent_ratio);
  return surface_integral_new(name, mlpg, vtable);
}

static int cube_vol_num_quad_points(void* context, int i)
{
  mlpg_simple_t* mlpg = context;
  return mlpg->N * mlpg->N * mlpg->N;
}

static void cube_vol_get_quad(void* context, int i, point_t* points, real_t* weights)
{
  mlpg_simple_t* mlpg = context;
  real_t gauss_pts[mlpg->N], gauss_wts[mlpg->N];
  get_gauss_legendre_points(mlpg->N, gauss_pts, gauss_wts);
  point_t* xi = &mlpg->cloud->points[i];
  real_t hi = mlpg->extents[i];
  real_t L = mlpg->ratio * hi;
  real_t a = xi->x - 0.5*L, b = xi->x + 0.5*L;
  real_t c = xi->y - 0.5*L, d = xi->y + 0.5*L;
  real_t e = xi->z - 0.5*L, f = xi->z + 0.5*L;
  real_t V = (b - a) * (d - c) * (f - e);
  int m = 0;
  for (int ii = 0; ii < mlpg->N; ++ii)
  {
    for (int jj = 0; jj < mlpg->N; ++jj)
    {
      for (int kk = 0; kk < mlpg->N; ++kk, ++m)
      {
        points[m].x = (b - a) * gauss_pts[ii] + a;
        points[m].y = (d - c) * gauss_pts[jj] + c;
        points[m].z = (f - e) * gauss_pts[kk] + e;
        weights[m] = gauss_wts[ii]*gauss_wts[jj]*gauss_wts[kk] * V;
      }
    }
  }
}

volume_integral_t* mlpg_cube_volume_integral_new(point_cloud_t* cloud,
                                                 real_t* extents,
                                                 int num_points,
                                                 real_t side_to_extent_ratio)
{
  mlpg_simple_t* mlpg = mlpg_simple_new(cloud, extents, num_points, 
                                        side_to_extent_ratio);
  volume_integral_vtable vtable = {.num_quad_points = cube_vol_num_quad_points,
{'fraction_non_alphanumeric': 0.15879065040650406, 'fraction_numerical': 0.02015582655826558, 'mean_word_length': 1.9858407079646017, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 91, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '325877', 'n_tokens_mistral': 5173, 'n_tokens_neox': 4997, 'n_words': 996}