Skip to content

GitLab

Commit 68db8d4e authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

mesh_find : minor_optimizations

parent 3bdc41c5
Hide whitespace changes
Inline Side-by-side
src/mesh_find.c
View file @ 68db8d4e
......@@ -4,10 +4,10 @@
#include "mesh.h"
#include "vector.h"
#include "mesh_find.h"
#define N (DIMENSION+1)
struct MeshTree_{
Cell *tree;
double *elements;
Mesh *mesh;
};
......@@ -23,7 +23,25 @@ static void _bboxtol(double *bbmin, double *bbmax) {
}
}
const double volTet(const double x0[3], const double x1[3], const double x2[3], const double x3[3]) {
#if DIMENSION == 2
static inline int isInsideEl(const double x[3], const double *X[3], double *Xi){
double dx = x[0] - X[0][0], dy = x[1] - X[0][1];
double DX[2] = {X[1][0] - X[0][0], X[2][0] - X[0][0]};
double DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
double det = DX[1] * DY[0] - DX[0] * DY[1];
double xi = (DX[1] * dy - DY[1] * dx) / det;
double eta = -(DX[0] * dy - DY[0] * dx) / det;
if (xi > -1e-8 && eta > -1e-8 && 1 - xi - eta > -1e-8) {
Xi[0] = xi;
Xi[1] = eta;
return 1;
}
return 0;
}
#else
static inline const double volTet(const double x0[3], const double x1[3], const double x2[3], const double x3[3]) {
// x1 * (y2 * (z3 - z4) - z2 * (y3 - y4) + y3 * z4 - y4 * z3)
//- y1 * (x2 * (z3 - z4) - z2 * (x3 - x4) + x3 * z4 - x4 * z3)
//+ z1 * (x2 * (y3 - y4) - y2 * (x3 - x4) + x3 * y4 - x4 * y3)
......@@ -33,24 +51,18 @@ const double volTet(const double x0[3], const double x1[3], const double x2[3],
+ x0[2] * (x1[0] * (x2[1] - x3[1]) - x1[1] * (x2[0] - x3[0]) + x2[0] * x3[1] - x3[0] * x2[1])
- (x1[0] * (x2[1] * x3[2] - x3[1] * x2[2]) - x1[1] * (x2[0] * x3[2] - x3[0] * x2[2]) + x1[2] * (x2[0] * x3[1] - x3[0] * x2[1]));
}
static inline int isInsideTet(const double x[3], const double x0[3], const double x1[3], const double x2[3], const double x3[3], double *xi){
for (int i = 0; i< 3; ++i) {
const double xx = x[i];
if ((xx < x0[i] && xx< x1[i] && xx < x2[i] && xx < x3[i])
||(xx > x0[i] && xx> x1[i] && xx > x2[i] && xx > x3[i]))
return 0;
}
double v = volTet(x0, x1, x2, x3);
double v0 = volTet(x, x1, x2, x3);
static inline int isInsideEl(const double x[3], const double *X[4], double *xi){
double v = volTet(X[0], X[1], X[2], X[3]);
double v0 = volTet(x, X[1], X[2], X[3]);
if (v0 / v < -1e-8)
return 0;
double v1 = volTet(x0, x, x2, x3);
double v1 = volTet(X[0], x, X[2], X[3]);
if (v1 / v < -1e-8)
return 0;
double v2 = volTet(x0, x1, x, x3);
double v2 = volTet(X[0], X[1], x, X[3]);
if (v2 / v < -1e-8)
return 0;
double v3 = volTet(x0, x1, x2, x);
double v3 = volTet(X[0], X[1], X[2], x);
if (v3 / v < -1e-8)
return 0;
......@@ -59,91 +71,39 @@ static inline int isInsideTet(const double x[3], const double x0[3], const doubl
xi[2] = v3 / v;
return 1;
}
#endif
void mesh_tree_particle_to_mesh(MeshTree *mt, size_t nParticles, double *particles, int *elid, double *Xi)
{
const Mesh *mesh = mt->mesh;
int *found = NULL;
vectorPushN(&found, 100);
vectorClear(found);
int N = DIMENSION + 1;
const double *elements = mt->elements;
for (size_t i = 0; i < nParticles; ++i) {
double *x = &particles[i * DIMENSION];
if (elid[i]!=-1) {
const double *el = &elements[elid[i] * N * DIMENSION];
#if DIMENSION == 2
const double *X[3] = {el, el + DIMENSION, el + DIMENSION * 2};
double dx = x[0] - X[0][0], dy = x[1] - X[0][1];
double DX[2] = {X[1][0] - X[0][0], X[2][0] - X[0][0]};
double DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
double det = DX[1] * DY[0] - DX[0] * DY[1];
double xi = (DX[1] * dy - DY[1] * dx) / det;
double eta = -(DX[0] * dy - DY[0] * dx) / det;
if (xi > -1e-8 && eta > -1e-8 && 1 - xi - eta > -1e-8) {
Xi[i * DIMENSION + 0] = xi;
Xi[i * DIMENSION + 1] = eta;
}
else
elid[i] = -1;
#else
const double *X[4] = {el, el + DIMENSION, el + DIMENSION * 2, el + DIMENSION * 3};
double xi[3];
if (!isInsideTet(x,el,el+3,el+6,el+9,Xi+i*DIMENSION))
const double *X[N];
for (int k = 0; k < N; ++k)
X[k] = mesh->x+mesh->elements[elid[i]*N+k]*3;
if (!isInsideEl(x,X,Xi+i*DIMENSION))
elid[i] = -1;
#endif
else
continue;
}
if (elid[i] != -1)
continue;
cellSearch(mt->tree, x, x, &found);
for (size_t j = 0; j < vectorSize(found); ++j) {
const double *el = &elements[found[j] * N * DIMENSION];
if (DIMENSION == 2) {
const double *X[3] = {el, el + DIMENSION, el + DIMENSION * 2};
double dx = x[0] - X[0][0], dy = x[1] - X[0][1];
double DX[2] = {X[1][0] - X[0][0], X[2][0] - X[0][0]};
double DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
double det = DX[1] * DY[0] - DX[0] * DY[1];
double xi = (DX[1] * dy - DY[1] * dx) / det;
double eta = -(DX[0] * dy - DY[0] * dx) / det;
if (xi > -1e-8 && eta > -1e-8 && 1 - xi - eta > -1e-8) {
Xi[i * DIMENSION + 0] = xi;
Xi[i * DIMENSION + 1] = eta;
elid[i] = found[j];
break;
}
}
else {
if (isInsideTet(x,el,el+3,el+6,el+9,Xi+i*DIMENSION)){
elid[i] = found[j];
break;
}
const double *X[N];
for (int k = 0; k < N; ++k)
X[k] = mesh->x+mesh->elements[found[j]*N+k]*3;
if (isInsideEl(x,X,Xi+i*DIMENSION)) {
elid[i] = found[j];
break;
}
}
/*if (elid[i] == -1) {
printf(" PARTICLE %i OUTSIDE DOMAIN N2 search!!!\n", i);
double toll = -1000;
for (size_t j = 0; j < nElements; ++j) {
double *el = &elements[j * N * DIMENSION];
double *X[3] = {el, el + DIMENSION, el + DIMENSION * 2};
double dx = x[0] - X[0][0], dy = x[1] - X[0][1];
double DX[2] = {X[1][0] - X[0][0], X[2][0] - X[0][0]};
double DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
double det = DX[1] * DY[0] - DX[0] * DY[1];
double xi = (DX[1] * dy - DY[1] * dx) / det;
double eta = -(DX[0] * dy - DY[0] * dx) / det;
toll = fmax(toll, fmin(-xi, fmin(-eta, xi + eta -1)));
if (xi > -1e-8 && eta > -1e-8 && 1 - xi - eta > -1e-8) {
Xi[i * DIMENSION + 0] = xi;
Xi[i * DIMENSION + 1] = eta;
elid[i] = j;
break;
}
}*/
if (elid[i] == -1) {
//printf(" PARTICLE %lu OUTSIDE DOMAIN!!! %g %g\n", i, x[0], x[1]);
//exit(1);
}
//}
if (elid[i] == -1) {
//printf(" PARTICLE %lu OUTSIDE DOMAIN!!! %g %g\n", i, x[0], x[1]);
//exit(1);
}
vectorClear(found);
}
vectorFree(found);
......@@ -153,49 +113,31 @@ void mesh_tree_particle_to_mesh(MeshTree *mt, size_t nParticles, double *particl
MeshTree *mesh_tree_create(Mesh *mesh) {
MeshTree *mt = malloc(sizeof(MeshTree));
mt->mesh = mesh;
size_t nElements = mesh->n_elements;
#if DIMENSION == 2
double *elements = malloc(sizeof(double)*6*mesh->n_elements);
for (int i = 0; i < nElements*3; ++i) {
int iv = mesh->elements[i];
elements[i*2+0] = mesh->x[iv*3+0];
elements[i*2+1] = mesh->x[iv*3+1];
}
#else
double *elements = malloc(sizeof(double)*12*mesh->n_elements);
for (int i = 0; i < mesh->n_elements*4; ++i) {
int iv = mesh->elements[i];
elements[i*3+0] = mesh->x[iv*3+0];
elements[i*3+1] = mesh->x[iv*3+1];
elements[i*3+2] = mesh->x[iv*3+2];
}
#endif
mt->elements = elements;
double bbmin[DIMENSION], bbmax[DIMENSION];
int N = DIMENSION + 1;
for (int i = 0; i < DIMENSION; ++i) {
bbmin[i] = elements[i];
bbmax[i] = elements[i];
bbmin[i] = mesh->x[i];
bbmax[i] = mesh->x[i];
}
for (int i = 0; i < N * nElements; ++i) {
for (int i = 0; i < mesh->n_nodes; ++i) {
for (int d = 0; d < DIMENSION; ++d) {
bbmin[d] = fmin(bbmin[d], elements[DIMENSION * i + d]);
bbmax[d] = fmax(bbmax[d], elements[DIMENSION * i + d]);
bbmin[d] = fmin(bbmin[d], mesh->x[3*i+d]);
bbmax[d] = fmax(bbmax[d], mesh->x[3*i+d]);
}
}
_bboxtol(bbmin, bbmax);
Cell *tree = cellNew(bbmin, bbmax);
double amin[DIMENSION], amax[DIMENSION];
for (int i = 0; i < nElements; ++i) {
double *el = elements + (DIMENSION * N) * i;
for (int i = 0; i < mesh->n_elements; ++i) {
const double *x = &mesh->x[mesh->elements[i*N+0]*3];
for (int d = 0; d < DIMENSION; ++d) {
amin[d] = el[d];
amax[d] = el[d];
amin[d] = x[d];
amax[d] = x[d];
}
for (int v = 1; v < N; ++v) {//only simplices
const double *x = &mesh->x[mesh->elements[i*N+v]*3];
for (int d = 0; d < DIMENSION; ++d) {
amin[d] = fmin(amin[d], el[v * DIMENSION + d]);
amax[d] = fmax(amax[d], el[v * DIMENSION + d]);
amin[d] = fmin(amin[d], x[d]);
amax[d] = fmax(amax[d], x[d]);
}
}
_bboxtol(amin, amax);
......@@ -205,9 +147,7 @@ MeshTree *mesh_tree_create(Mesh *mesh) {
return mt;
}
void mesh_tree_free(MeshTree *mt) {
cellFree(mt->tree);
free(mt->elements);
free(mt);
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment