Skip to content
Snippets Groups Projects
Commit 04da0ca3 authored by Hans Stephensen's avatar Hans Stephensen
Browse files

added c++ code and swig interface for python

parent 6a7d1ec3
Branches master
No related tags found
No related merge requests found
Showing
with 7491 additions and 1 deletion
%module NeighbourhoodTetMeshIntegrator
%{
#define SWIG_FILE_WITH_INIT
#include "c_NeighbourhoodTetMeshIntegrator.h"
%}
%include "numpy.i"
%init %{
import_array();
%}
%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) {
(double * verts, int num_verts, int d1)
};
%apply (long* INPLACE_ARRAY2, int DIM1, int DIM2) {
(const long * const elements, int num_elements, int d2)
};
%apply (double* INPLACE_ARRAY1, int DIM1) {
(double * dists, int num_dists),
(const double * const x, int num_x),
(double * y_volume, int num_y1),
(double * y_area, int num_y2)
};
void NeighbourhoodTetMeshIntegrator(
double * verts, int num_verts, int d1,
const long * const elements, int num_elements, int d2,
double * dists, int num_dists,
const double * const x, int num_x,
double * y_volume, int num_y1,
double * y_area, int num_y2
);
\ No newline at end of file
# This file was automatically generated by SWIG (http://www.swig.org).
# Version 3.0.12
#
# Do not make changes to this file unless you know what you are doing--modify
# the SWIG interface file instead.
from sys import version_info as _swig_python_version_info
if _swig_python_version_info >= (2, 7, 0):
def swig_import_helper():
import importlib
pkg = __name__.rpartition('.')[0]
mname = '.'.join((pkg, '_NeighbourhoodTetMeshIntegrator')).lstrip('.')
try:
return importlib.import_module(mname)
except ImportError:
return importlib.import_module('_NeighbourhoodTetMeshIntegrator')
_NeighbourhoodTetMeshIntegrator = swig_import_helper()
del swig_import_helper
elif _swig_python_version_info >= (2, 6, 0):
def swig_import_helper():
from os.path import dirname
import imp
fp = None
try:
fp, pathname, description = imp.find_module('_NeighbourhoodTetMeshIntegrator', [dirname(__file__)])
except ImportError:
import _NeighbourhoodTetMeshIntegrator
return _NeighbourhoodTetMeshIntegrator
try:
_mod = imp.load_module('_NeighbourhoodTetMeshIntegrator', fp, pathname, description)
finally:
if fp is not None:
fp.close()
return _mod
_NeighbourhoodTetMeshIntegrator = swig_import_helper()
del swig_import_helper
else:
import _NeighbourhoodTetMeshIntegrator
del _swig_python_version_info
try:
_swig_property = property
except NameError:
pass # Python < 2.2 doesn't have 'property'.
try:
import builtins as __builtin__
except ImportError:
import __builtin__
def _swig_setattr_nondynamic(self, class_type, name, value, static=1):
if (name == "thisown"):
return self.this.own(value)
if (name == "this"):
if type(value).__name__ == 'SwigPyObject':
self.__dict__[name] = value
return
method = class_type.__swig_setmethods__.get(name, None)
if method:
return method(self, value)
if (not static):
if _newclass:
object.__setattr__(self, name, value)
else:
self.__dict__[name] = value
else:
raise AttributeError("You cannot add attributes to %s" % self)
def _swig_setattr(self, class_type, name, value):
return _swig_setattr_nondynamic(self, class_type, name, value, 0)
def _swig_getattr(self, class_type, name):
if (name == "thisown"):
return self.this.own()
method = class_type.__swig_getmethods__.get(name, None)
if method:
return method(self)
raise AttributeError("'%s' object has no attribute '%s'" % (class_type.__name__, name))
def _swig_repr(self):
try:
strthis = "proxy of " + self.this.__repr__()
except __builtin__.Exception:
strthis = ""
return "<%s.%s; %s >" % (self.__class__.__module__, self.__class__.__name__, strthis,)
try:
_object = object
_newclass = 1
except __builtin__.Exception:
class _object:
pass
_newclass = 0
def NeighbourhoodTetMeshIntegrator(verts, elements, dists, x, y_volume, y_area):
return _NeighbourhoodTetMeshIntegrator.NeighbourhoodTetMeshIntegrator(verts, elements, dists, x, y_volume, y_area)
NeighbourhoodTetMeshIntegrator = _NeighbourhoodTetMeshIntegrator.NeighbourhoodTetMeshIntegrator
# This file is compatible with both classic and new-style classes.
import numpy as np
from .NeighbourhoodTetMeshIntegrator import NeighbourhoodTetMeshIntegrator
def TetMeshIntegrator(verts, elements, dists, x):
verts = verts.astype(np.float64)
elements = elements.astype(np.int64)
dists = dists.astype(np.float64)
x = x.astype(np.float64)
y_volume = np.empty(x.shape, dtype=np.float64)
y_area = np.empty(x.shape, dtype=np.float64)
NeighbourhoodTetMeshIntegrator(verts, elements, dists, x, y_volume, y_area)
return y_volume, y_area
\ No newline at end of file
File added
from .TetMeshIntegrator import TetMeshIntegrator
\ No newline at end of file
#include "c_NeighbourhoodTetMeshIntegrator.h"
#include <cmath>
#include <algorithm>
#include <iostream>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
#define TOL 1e-10
void copy_verts(const double * const verts, const long * const v_indices, double * t_verts, const int * const ordering)
{
for (int j = 0; j < 4; ++j)
{
long vn = v_indices[ordering[j]];
for (int i = 0; i < 3; ++i)
{
t_verts[3*j + i] = verts[3*vn + i];
}
}
}
void get_sorted_indices(double * dists, int * ordering)
{
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3-i; ++j) // not sure about this, verify with web for bubblesort
{
if (dists[j] > dists[j+1]) {
std::swap(dists[j], dists[j+1]);
std::swap(ordering[j], ordering[j+1]);
}
}
}
void rotx(double & vy, double & vz, const double t)
{
const double ct = std::cos(t);
const double st = std::sin(t);
const double tmp_vz = vy * st + vz * ct;
vy = vy * ct - vz * st;
vz = tmp_vz;
}
void roty(double & vx, double & vz, const double t)
{
const double ct = std::cos(t);
const double st = std::sin(t);
const double tmp_vz = ct * vz - st * vx;
vx = vx * ct + st * vz;
vz = tmp_vz;
}
void rotz(double & vx, double & vy, const double t)
{
const double ct = std::cos(t);
const double st = std::sin(t);
const double tmp_vy = vx * st + vy * ct;
vx = vx * ct - vy * st;
vy = tmp_vy;
}
void get_linear_params(const double * const verts, const double * const dists, double &a, double &b, double &c)
{
using namespace boost::numeric::ublas;
matrix<double> A(4,4);
vector<double> v(4);
for (unsigned i = 0; i < 4; ++i) {
v(i) = dists[i];
for (unsigned j = 0; j < 4; ++j) {
if (j < 3) A(i, j) = verts[i*3 + j];
else A(i, j) = 1.0;
}
}
permutation_matrix<size_t> pm(4);
int singular = lu_factorize(A, pm);
if (singular) {
a=0; b=0; c=0; // This is not my favorite solution, but it is caught later down the line
return;
}
lu_substitute(A, pm, v);
a = v(0);
b = v(1);
c = v(2);
}
double full_tri_area(const double * const verts)
{
// cross product and multiply by last vector
double a[3] = {verts[3] - verts[0], verts[4] - verts[1], verts[5] - verts[2]};
double b[3] = {verts[6] - verts[0], verts[7] - verts[1], verts[8] - verts[2]};
double cx = a[1] * b[2] - b[1] * a[2];
double cy = b[0] * a[2] - a[0] * b[2];
double cz = a[0] * b[1] - b[0] * a[1];
return std::sqrt(cx*cx + cy*cy + cz*cz) / 2.0;
}
double full_tet_volume(const double * const verts)
{
// cross product and multiply by last vector
double a[3] = {verts[3] - verts[0], verts[4] - verts[1], verts[5] - verts[2]};
double b[3] = {verts[6] - verts[0], verts[7] - verts[1], verts[8] - verts[2]};
double c[3] = {verts[9] - verts[0], verts[10] - verts[1], verts[11] - verts[2]};
double cm_x = (b[1] * c[2] - c[1] * b[2]) * a[0];
double cm_y = (c[0] * b[2] - b[0] * c[2]) * a[1];
double cm_z = (b[0] * c[1] - c[0] * b[1]) * a[2];
return std::abs(cm_x + cm_y + cm_z) / 6.0;
}
void transform_tet(double * verts, double * dists, double & scale) // 9, 3
{
// Normalize so first vertice is at the origin
for (int j = 1; j < 4; ++j)
for (int i = 0; i < 3; ++i)
verts[j*3 + i] -= verts[i];
for (int i = 0; i < 3; ++i)
verts[i] = 0;
double a, b, c;
get_linear_params(verts, dists, a, b, c);
double gradient_norm = std::sqrt(a*a + b*b + c*c);
scale = gradient_norm;
double t;
t = -std::atan2(c, b);
for (int i = 1; i < 4; ++i)
{
// rotx (x is unchanged, rot y and z)
rotx(verts[i*3+1], verts[i*3+2], t);
}
rotx(b, c, t);
t = -std::atan2(b, a);
for (int i = 1; i < 4; ++i)
{
// rotz (z is unchanged, rot x and y)
rotz(verts[i*3], verts[i*3+1], t);
}
// Scale to h scale
for (int i = 0; i < 4; ++i)
{
verts[i*3] = verts[i*3] * gradient_norm + dists[0];
}
}
void h_midpoint(const double * const p0, const double * const p1, const double h, double * p01)
{
double K[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]};
double h0 = p0[0];
double h1 = p1[0];
double a = h0;
double b = 1.0/(h1 - h0);
for (int i = 0; i < 3; ++i)
{
p01[i] = p0[i] + K[i]*(h - a)*b;
}
}
double volume_fast(const double * const v0, const double * const v1,
const double * const v2, const double * const v3,
const double h)
{
double p01[3];
h_midpoint(v0, v1, h, p01);
double p02[3];
h_midpoint(v0, v2, h, p02);
double p03[3];
h_midpoint(v0, v3, h, p03);
double p12[3];
h_midpoint(v1, v2, h, p12);
double p13[3];
h_midpoint(v1, v3, h, p13);
double P_tet1[12] = { v0[0], v0[1], v0[2],
p01[0], p01[1], p01[2],
p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2]};
double P_tet2[12] = { v1[0], v1[1], v1[2],
p01[0], p01[1], p01[2],
p12[0], p12[1], p12[2],
p13[0], p13[1], p13[2]};
double volume_1 = full_tet_volume(P_tet1);
double volume_2 = full_tet_volume(P_tet2);
return volume_1 - volume_2;
}
double volume_safe(const double * const v0, const double * const v1,
const double * const v2, const double * const v3,
const double h)
{
double p02[3];
h_midpoint(v0, v2, h, p02);
double p03[3];
h_midpoint(v0, v3, h, p03);
double p12[3];
h_midpoint(v1, v2, h, p12);
double p13[3];
h_midpoint(v1, v3, h, p13);
double P_tet1[12] = { v0[0], v0[1], v0[2],
p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2],
p12[0], p12[1], p12[2]};
double P_tet2[12] = { v0[0], v0[1], v0[2],
p03[0], p03[1], p03[2],
p12[0], p12[1], p12[2],
p13[0], p13[1], p13[2]};
double P_tet3[12] = { v0[0], v0[1], v0[2],
v1[0], v1[1], v1[2],
p12[0], p12[1], p12[2],
p13[0], p13[1], p13[2]};
double volume_1 = full_tet_volume(P_tet1);
double volume_2 = full_tet_volume(P_tet2);
double volume_3 = full_tet_volume(P_tet3);
return volume_1 + volume_2 + volume_3;
}
double sliced_tet_volume(const double * const verts, const double full_volume, const int i, double h)
{
double v0[3] = {verts[i*12 + 0], verts[i*12 + 1], verts[i*12 + 2]};
double v1[3] = {verts[i*12 + 3], verts[i*12 + 4], verts[i*12 + 5]};
double v2[3] = {verts[i*12 + 6], verts[i*12 + 7], verts[i*12 + 8]};
double v3[3] = {verts[i*12 + 9], verts[i*12 + 10], verts[i*12 + 11]};
if (h <= v0[0]) {
// VOL 0
return 0.0;
} else if (h <= v1[0]) {
// VOL 1
double p01[3];
double p02[3];
double p03[3];
h_midpoint(v0, v1, h, p01);
h_midpoint(v0, v2, h, p02);
h_midpoint(v0, v3, h, p03);
double P[12] = { v0[0], v0[1], v0[2],
p01[0], p01[1], p01[2],
p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2]};
double volume = full_tet_volume(P);
if (std::isnan(volume))
{
std::cout << "NAN in volume (case 1) on " << P << " " << h << std::endl;
}
return volume;
} else if (h <= v2[0]) {
// VOL 2
if (std::abs(v0[0] - v1[0]) < TOL) {
// Degenerate case
double volume = volume_safe(v0, v1, v2, v3, h);
if (std::isnan(volume))
{
std::cout << "NAN in volume_safe on " << v0 << " " << v1 << " " << v2 << " " << v3 << " " << h << std::endl;
}
return volume;
} else {
// Non-degenerate case
double volume = volume_fast(v0, v1, v2, v3, h);
if (std::isnan(volume))
{
std::cout << "NAN in volume_fast on " << v0 << " " << v1 << " " << v2 << " " << v3 << " " << h << std::endl;
}
return volume;
}
} else if (h <= v3[0]) {
// VOL 3
double p03[3];
double p13[3];
double p23[3];
h_midpoint(v0, v3, h, p03);
h_midpoint(v1, v3, h, p13);
h_midpoint(v2, v3, h, p23);
double P_remaining[12] = { v3[0], v3[1], v3[2],
p03[0], p03[1], p03[2],
p13[0], p13[1], p13[2],
p23[0], p23[1], p23[2]};
double volume_remaining = full_tet_volume(P_remaining);
if (std::isnan(volume_remaining))
{
std::cout << "NAN in volume_remaining on " << P_remaining << " " << h << std::endl;
}
if (std::isnan(full_volume))
{
std::cout << "NAN in full_volume on " << v0 << " " << v1 << " " << v2 << " " << v3 << " " << h << std::endl;
}
return full_volume - volume_remaining;
} else {
// VOL 4
return full_volume;
}
}
double sliced_tet_area(const double * const verts, const int i, double h)
{
double v0[3] = {verts[i*12 + 0], verts[i*12 + 1], verts[i*12 + 2]};
double v1[3] = {verts[i*12 + 3], verts[i*12 + 4], verts[i*12 + 5]};
double v2[3] = {verts[i*12 + 6], verts[i*12 + 7], verts[i*12 + 8]};
double v3[3] = {verts[i*12 + 9], verts[i*12 + 10], verts[i*12 + 11]};
if (h <= v0[0]) {
// AREA 0
return 0.0;
} else if (h <= v1[0]) {
// AREA 1
double p01[3];
double p02[3];
double p03[3];
h_midpoint(v0, v1, h, p01);
h_midpoint(v0, v2, h, p02);
h_midpoint(v0, v3, h, p03);
double P[9] = {p01[0], p01[1], p01[2],
p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2]};
return full_tri_area(P);
} else if (h <= v2[0]) {
// AREA 2
double p02[3];
double p03[3];
double p12[3];
double p13[3];
h_midpoint(v0, v2, h, p02);
h_midpoint(v0, v3, h, p03);
h_midpoint(v1, v2, h, p12);
h_midpoint(v1, v3, h, p13);
double P_1[9] = {p02[0], p02[1], p02[2],
p03[0], p03[1], p03[2],
p12[0], p12[1], p12[2]};
double P_2[9] = {p03[0], p03[1], p03[2],
p12[0], p12[1], p12[2],
p13[0], p13[1], p13[2]};
return full_tri_area(P_1) + full_tri_area(P_2);
} else if (h <= v3[0]) {
// AREA 3
double p03[3];
double p13[3];
double p23[3];
h_midpoint(v0, v3, h, p03);
h_midpoint(v1, v3, h, p13);
h_midpoint(v2, v3, h, p23);
double P[9] = {p03[0], p03[1], p03[2],
p13[0], p13[1], p13[2],
p23[0], p23[1], p23[2]};
return full_tri_area(P);
} else {
// AREA 4
return 0.0;
}
}
void NeighbourhoodTetMeshIntegrator(
double * verts, int num_verts, int d1,
const long * const elements, int num_elements, int d2,
double * dists, int num_dists,
const double * const x, int num_x,
double * y_volume, int num_y1,
double * y_area, int num_y2
)
{
double * t_verts = (double *) malloc(num_elements * 12 * sizeof(double)); // Every triangle, x0 y0 x1 y1 x2 y2
double * t_scales = (double *) malloc(num_elements * sizeof(double));
double * full_volumes = (double *) malloc(num_elements * sizeof(double));
#pragma omp parallel for
for (int i = 0; i < num_elements; ++i)
{
double * t_verts_tmp = (double *) malloc(12*sizeof(double));
int ordering[4];
long v_indices[4];
double t_dists_tmp[4];
for (int j = 0; j < 4; ++j)
{
ordering[j] = j;
v_indices[j] = elements[i*4+j];
t_dists_tmp[j] = dists[v_indices[j]];
}
get_sorted_indices(t_dists_tmp, ordering);
copy_verts(verts, v_indices, t_verts_tmp, ordering);
double t_scale;
transform_tet(t_verts_tmp, t_dists_tmp, t_scale);
t_scales[i] = t_scale;
full_volumes[i] = full_tet_volume(t_verts_tmp);
for (int j = 0; j < 12; ++j)
{
t_verts[i*12 + j] = t_verts_tmp[j];
}
free(t_verts_tmp);
}
#pragma omp parallel for
for (int j = 0; j < num_x; ++j)
{
y_volume[j] = 0.0;
y_area[j] = 0.0;
for (int i = 0; i < num_elements; ++i)
{
if (t_scales[i] != 0)
{
y_volume[j] += sliced_tet_volume(t_verts, full_volumes[i], i, x[j]) / t_scales[i];
y_area[j] += sliced_tet_area(t_verts, i, x[j]);
}
}
}
free(t_verts);
free(t_scales);
free(full_volumes);
}
\ No newline at end of file
#ifndef C_NEIGHBOURHOODTETMESHINTEGRATOR_HPP
#define C_NEIGHBOURHOODTETMESHINTEGRATOR_HPP
#define PI 3.141592653589793238462643383279
void NeighbourhoodTetMeshIntegrator(
double * verts, int num_verts, int d1,
const long * const elements, int num_elements, int d2,
double * dists, int num_dists,
const double * const x, int num_x,
double * y_volume, int num_y1,
double * y_area, int num_y2
);
#endif // C_NEIGHBOURHOODTETMESHINTEGRATOR_HPP
\ No newline at end of file
This diff is collapsed.
swig -Wall -python -c++ NeighbourhoodTetMeshIntegrator.i
g++ -O3 -fPIC -fopenmp -lstdc++ -c c_NeighbourhoodTetMeshIntegrator.cpp
g++ -O3 -fPIC -fopenmp -lstdc++ -c NeighbourhoodTetMeshIntegrator_wrap.cxx -I/usr/include/python3.7 -I/usr/lib/python3.7/site-packages/numpy/core/include
g++ -shared c_NeighbourhoodTetMeshIntegrator.o NeighbourhoodTetMeshIntegrator_wrap.o -fopenmp -O3 -lstdc++ -o _NeighbourhoodTetMeshIntegrator.so
%module NeighbourhoodTriMeshIntegrator
%{
#define SWIG_FILE_WITH_INIT
#include "c_NeighbourhoodTriMeshIntegrator.h"
%}
%include "numpy.i"
%init %{
import_array();
%}
%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) {
(double * verts, int num_verts, int d1)
};
%apply (long* INPLACE_ARRAY2, int DIM1, int DIM2) {
(long * faces, int num_faces, int d2)
};
%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) {
(double * tet_verts, int num_tet_verts, int d3)
};
%apply (long* INPLACE_ARRAY2, int DIM1, int DIM2) {
(long * tet_elems, int num_tet_verts, int d4)
};
%apply (long* INPLACE_ARRAY2, int DIM1, int DIM2) {
(long * tet_faces, int num_tet_verts, int d5)
};
%apply (double* INPLACE_ARRAY1, int DIM1) {
(double * dists, int num_dists),
(double * x, int num_x),
(double * y_area, int num_y1),
(double * y_circumference, int num_y2)
};
void NeighbourhoodTriMeshIntegrator(
double * verts, int num_verts, int d1,
long * faces, int num_faces, int d2,
double * dists, int num_dists,
double * x, int num_x,
double * y_area, int num_y1,
double * y_circumference, int num_y2
);
# This file was automatically generated by SWIG (http://www.swig.org).
# Version 4.0.2
#
# Do not make changes to this file unless you know what you are doing--modify
# the SWIG interface file instead.
from sys import version_info as _swig_python_version_info
if _swig_python_version_info < (2, 7, 0):
raise RuntimeError("Python 2.7 or later required")
# Import the low-level C/C++ module
if __package__ or "." in __name__:
from . import _NeighbourhoodTriMeshIntegrator
else:
import _NeighbourhoodTriMeshIntegrator
try:
import builtins as __builtin__
except ImportError:
import __builtin__
def _swig_repr(self):
try:
strthis = "proxy of " + self.this.__repr__()
except __builtin__.Exception:
strthis = ""
return "<%s.%s; %s >" % (self.__class__.__module__, self.__class__.__name__, strthis,)
def _swig_setattr_nondynamic_instance_variable(set):
def set_instance_attr(self, name, value):
if name == "thisown":
self.this.own(value)
elif name == "this":
set(self, name, value)
elif hasattr(self, name) and isinstance(getattr(type(self), name), property):
set(self, name, value)
else:
raise AttributeError("You cannot add instance attributes to %s" % self)
return set_instance_attr
def _swig_setattr_nondynamic_class_variable(set):
def set_class_attr(cls, name, value):
if hasattr(cls, name) and not isinstance(getattr(cls, name), property):
set(cls, name, value)
else:
raise AttributeError("You cannot add class attributes to %s" % cls)
return set_class_attr
def _swig_add_metaclass(metaclass):
"""Class decorator for adding a metaclass to a SWIG wrapped class - a slimmed down version of six.add_metaclass"""
def wrapper(cls):
return metaclass(cls.__name__, cls.__bases__, cls.__dict__.copy())
return wrapper
class _SwigNonDynamicMeta(type):
"""Meta class to enforce nondynamic attributes (no new attributes) for a class"""
__setattr__ = _swig_setattr_nondynamic_class_variable(type.__setattr__)
def NeighbourhoodTriMeshIntegrator(verts, faces, dists, x, y_area, y_circumference):
return _NeighbourhoodTriMeshIntegrator.NeighbourhoodTriMeshIntegrator(verts, faces, dists, x, y_area, y_circumference)
import numpy as np
from .NeighbourhoodTriMeshIntegrator import NeighbourhoodTriMeshIntegrator
def TriMeshIntegrator(verts, faces, dists, x):
verts = verts.astype(np.float64)
faces = faces.astype(np.int64)
dists = dists.astype(np.float64)
x = x.astype(np.float64)
y_area = np.empty(x.shape, dtype=np.float64)
y_circumference = np.empty(x.shape, dtype=np.float64)
NeighbourhoodTriMeshIntegrator(verts, faces, dists, x, y_area, y_circumference)
return y_area, y_circumference
\ No newline at end of file
File added
from .TriMeshIntegrator import TriMeshIntegrator
\ No newline at end of file
#include "c_NeighbourhoodTriMeshIntegrator.h"
#include <cmath>
#include <algorithm>
#include <iostream>
#include <vector>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
const double eps = 0.000000001;
void copy_verts(const double * const verts, const long * const v_indices, double * t_verts, const int * const ordering)
{
for (int j = 0; j < 3; ++j)
{
long vn = v_indices[ordering[j]];
for (int i = 0; i < 3; ++i)
{
t_verts[3*j + i] = verts[3*vn + i];
}
}
}
void get_sorted_indices(double * dists, int * ordering)
{
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 2-i; ++j)
{
if (dists[j] > dists[j+1]) {
std::swap(dists[j], dists[j+1]);
std::swap(ordering[j], ordering[j+1]);
}
}
}
void rotx(double & vy, double & vz, const double t)
{
const double ct = std::cos(t);
const double st = std::sin(t);
const double tmp_vz = vy * st + vz * ct;
vy = vy * ct - vz * st;
vz = tmp_vz;
}
void roty(double & vx, double & vz, const double t)
{
const double ct = std::cos(t);
const double st = std::sin(t);
const double tmp_vz = ct * vz - st * vx;
vx = vx * ct + st * vz;
vz = tmp_vz;
}
void rotz(double & vx, double & vy, const double t)
{
const double ct = std::cos(t);
const double st = std::sin(t);
const double tmp_vy = vx * st + vy * ct;
vx = vx * ct - vy * st;
vy = tmp_vy;
}
void get_linear_params(const double * const verts, const double * const dists, double &a, double &b)
{
using namespace boost::numeric::ublas;
matrix<double> A(3,3);
vector<double> v(3);
for (unsigned i = 0; i < 3; ++i) {
v(i) = dists[i];
for (unsigned j = 0; j < 3; ++j) {
if (j < 2) A(i, j) = verts[i*3 + j];
else A(i, j) = 1.0;
}
}
permutation_matrix<size_t> pm(3);
lu_factorize(A, pm);
lu_substitute(A, pm, v);
a = v(0);
b = v(1);
}
void get_edge_linear_params(const double * const verts, double * linear_params, int i)
{
const double x0 = verts[i*6];
const double y0 = verts[i*6 + 1];
const double x1 = verts[i*6 + 2];
const double y1 = verts[i*6 + 3];
const double x2 = verts[i*6 + 4];
const double y2 = verts[i*6 + 5];
linear_params[i*6] = (y1 - y0) / (x1 - x0);
linear_params[i*6 + 1] = y0 - linear_params[i*6]*x0;
linear_params[i*6 + 2] = (y2 - y0) / (x2 - x0);
linear_params[i*6 + 3] = y0 - linear_params[i*6+2]*x0;
linear_params[i*6 + 4] = (y2 - y1) / (x2 - x1);
linear_params[i*6 + 5] = y1 - linear_params[i*6+4]*x1;
}
void transform_triangle(double * verts, double * dists, double & scale, double & full_area) // 9, 3
{
// Normalize so first vertice is at the origin
for (int i = 0; i < 3; ++i)
{
verts[3+i] -= verts[i];
verts[6+i] -= verts[i];
verts[i] = 0;
}
// Rotate so it corresponds to a flat triangWle
double t;
t = -std::atan2(verts[4], verts[3]);
rotz(verts[3],verts[4], t);
rotz(verts[6],verts[7], t);
t = std::atan2(verts[5], verts[3]);
roty(verts[3],verts[5], t);
roty(verts[6],verts[8], t);
t = -std::atan2(verts[8], verts[7]);
rotx(verts[4],verts[5], t);
rotx(verts[7],verts[8], t);
// calculate full area of triangle
const double ax = verts[0] - verts[3];
const double ay = verts[1] - verts[4];
const double bx = verts[0] - verts[6];
const double by = verts[1] - verts[7];
full_area = std::abs(ax*by - bx*ay) / 2.0;
double a, b;
if (full_area < eps) {
verts[0] = dists[0]; verts[1] = dists[0]; verts[2] = dists[0]; // dont think we need this now we have the bounds elsewhere (Update: well, apparently we do...)
return;
}
if (dists[0] == dists[2]) {
verts[0] = dists[0]; verts[1] = dists[0]; verts[2] = dists[0]; // dont think we need this now we have the bounds elsewhere (Update: well, apparently we do...)
return;
}
get_linear_params(verts, dists, a, b);
double gradient_norm = std::sqrt(a*a + b*b);
scale = 1.0/gradient_norm;
// Change basis
t = -std::atan2(b, a);
rotz(verts[3],verts[4], t);
rotz(verts[6],verts[7], t);
// Scale to h scale
verts[0] = verts[0] + dists[0];
verts[3] = verts[3] * gradient_norm + dists[0];
verts[6] = verts[6] * gradient_norm + dists[0];
}
double sliced_triangle_area(const double * const verts, const double * const edge_linear_params, const int i, double h)
{
double x0 = verts[i*6 + 0];
double y0 = verts[i*6 + 1];
double x1 = verts[i*6 + 2];
double y1 = verts[i*6 + 3];
double x2 = verts[i*6 + 4];
double y2 = verts[i*6 + 5];
auto f01 = [=](double x) { return edge_linear_params[i*6 + 0]*x + edge_linear_params[i*6 + 1]; };
auto f02 = [=](double x) { return edge_linear_params[i*6 + 2]*x + edge_linear_params[i*6 + 3]; };
auto f12 = [=](double x) { return edge_linear_params[i*6 + 4]*x + edge_linear_params[i*6 + 5]; };
if (x1 == x2 && y1 == y2) return 0.0;
if (x0 == x1) {
if (y0 == y1) return 0.0;
if (x0 == x2) return 0.0;
double height_full = x2 - x0;
double base_full = y1 - y0;
double full_area = std::abs(0.5*height_full*base_full);
double right_height = x2 - h;
double right_base = f02(h) - f12(h);
double right_area = std::abs(0.5*right_height*right_base);
double area = full_area - right_area;
return area;
}
if (h < x1)
{
double height = h - x0;
double base = f01(h) - f02(h);
double area = std::abs(0.5*base*height);
return std::abs(0.5*base*height);
}
double left_height = x1 - x0;
double base = f01(x1) - f02(x1);
double left_area = std::abs(base*left_height);
double right_height = x2 - x1;
double mid_base = y1 - f02(x1);
double full_right_area = std::abs(right_height*mid_base);
double remaining_right_area = std::abs((f12(h) - f02(h)) * (x2 - h));
double area = 0.5*(left_area + full_right_area - remaining_right_area);
return area;
}
double sliced_triangle_circ(const double * const verts, const double * const edge_linear_params, const int i, double h)
{
double x0 = verts[i*6 + 0];
double y0 = verts[i*6 + 1];
double x1 = verts[i*6 + 2];
double y1 = verts[i*6 + 3];
double x2 = verts[i*6 + 4];
double y2 = verts[i*6 + 5];
if (x0 == x1 && y0 == y1) return 0.0;
if (x0 == x2 && y0 == y2) return 0.0;
if (x1 == x2 && y1 == y2) return 0.0;
if (x0 == x1 && x1 == x2) return 0.0;
if (y0 == y1 && y1 == y2) return 0.0;
auto f02 = [=](double x) { return edge_linear_params[i*6 + 2]*x + edge_linear_params[i*6 + 3]; };
if (h < x1)
{
auto f01 = [=](double x) { return edge_linear_params[i*6 + 0]*x + edge_linear_params[i*6 + 1]; };
return std::abs(f01(h) - f02(h));
}
auto f12 = [=](double x) { return edge_linear_params[i*6 + 4]*x + edge_linear_params[i*6 + 5]; };
return std::abs(f12(h) - f02(h));
}
void NeighbourhoodTriMeshIntegrator(
double * verts, int num_verts, int d1,
long * faces, int num_faces, int d2,
double * dists, int num_dists,
double * x, int num_x,
double * y_area, int num_y1,
double * y_circumference, int num_y2
)
{
double * t_verts = (double *) malloc(num_faces * 6 * sizeof(double)); // Every triangle, x0 y0 x1 y1 x2 y2
double * t_scales = (double *) malloc(num_faces * sizeof(double));
double * full_areas = (double *) malloc(num_faces * sizeof(double));
double * triangle_h_bounds = (double *) malloc(num_faces * 2 * sizeof(double));
double * edge_linear_params = (double *) malloc (num_faces * 6 * sizeof(double));
#pragma omp parallel for
for (int i = 0; i < num_faces; ++i)
{
double * t_verts_tmp = (double *) malloc(9*sizeof(double));
int ordering[3] = {0, 1, 2};
long v_indices[3] = {faces[i*3 + 0], faces[i*3 + 1], faces[i*3 + 2]};
double t_dists_tmp[3] = {dists[v_indices[0]], dists[v_indices[1]], dists[v_indices[2]]};
get_sorted_indices(t_dists_tmp, ordering);
copy_verts(verts, v_indices, t_verts_tmp, ordering);
triangle_h_bounds[i*2] = t_dists_tmp[0];
triangle_h_bounds[i*2+1] = t_dists_tmp[2];
double t_scale, full_area;
transform_triangle(t_verts_tmp, t_dists_tmp, t_scale, full_area);
for (int j = 0; j < 6; ++j)
{
int g_idx = (j / 2) * 3 + j%2; // This skips the z-coordinate which is now 0
t_verts[i*6 + j] = t_verts_tmp[g_idx];
}
t_scales[i] = t_scale;
full_areas[i] = full_area;
get_edge_linear_params(t_verts, edge_linear_params, i);
free(t_verts_tmp);
}
#pragma omp parallel for
for (int j = 0; j < num_x; ++j)
{
y_area[j] = 0.0;
y_circumference[j] = 0.0;
for (int i = 0; i < num_faces; ++i)
{
if (x[j] <= triangle_h_bounds[i*2]) {
} else if (x[j] >= triangle_h_bounds[i*2+1]) {
y_area[j] += full_areas[i];
} else {
y_area[j] += sliced_triangle_area(t_verts, edge_linear_params, i, x[j]) * t_scales[i];
y_circumference[j] += sliced_triangle_circ(t_verts, edge_linear_params, i, x[j]);
}
}
}
free(t_verts);
free(t_scales);
free(full_areas);
free(triangle_h_bounds);
free(edge_linear_params);
}
#ifndef C_NEIGHBOURHOODTRIMESHINTEGRATOR_HPP
#define C_NEIGHBOURHOODTRIMESHINTEGRATOR_HPP
//#define PI 3.141592653589793238462643383279
void NeighbourhoodTriMeshIntegrator(
double * verts, int num_verts, int d1,
long * faces, int num_faces, int d2,
double * dists, int num_dists,
double * x, int num_x,
double * y_area, int num_y1,
double * y_circumference, int num_y2
);
#endif // C_NEIGHBOURHOODTRIMESHINTEGRATOR_HPP
\ No newline at end of file
This diff is collapsed.
swig -Wall -python -c++ NeighbourhoodTriMeshIntegrator.i
g++ -O3 -fPIC -fopenmp -lstdc++ -c c_NeighbourhoodTriMeshIntegrator.cpp
g++ -O3 -fPIC -fopenmp -lstdc++ -c NeighbourhoodTriMeshIntegrator_wrap.cxx -I/usr/include/python3.7 -I/usr/lib/python3.7/site-packages/numpy/core/include
g++ -shared c_NeighbourhoodTriMeshIntegrator.o NeighbourhoodTriMeshIntegrator_wrap.o -fopenmp -O3 -lstdc++ -o _NeighbourhoodTriMeshIntegrator.so
...@@ -2,7 +2,7 @@ import numpy as np ...@@ -2,7 +2,7 @@ import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from RandomRotations import RandomRotationMatrix from RandomRotations import RandomRotationMatrix
from RPSM import measure_interior, measure_exterior from src import measure_interior, measure_exterior
def Test_Cube_Integration_Interior(do_plots): def Test_Cube_Integration_Interior(do_plots):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment