Commit 3bdc41c5 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

fix bug : mesh to mesh projection much faster

parent c62e7f3b
...@@ -14,6 +14,7 @@ typedef struct { ...@@ -14,6 +14,7 @@ typedef struct {
int particle[MAXP]; int particle[MAXP];
double xmin[MAXP][DIMENSION]; double xmin[MAXP][DIMENSION];
double xmax[MAXP][DIMENSION]; double xmax[MAXP][DIMENSION];
int n;
} LeafData; } LeafData;
struct Cell_{ struct Cell_{
...@@ -31,8 +32,7 @@ void cellInit(Cell *c, const double *xmin, const double *xmax) { ...@@ -31,8 +32,7 @@ void cellInit(Cell *c, const double *xmin, const double *xmax) {
} }
c->sub = NULL; c->sub = NULL;
c->leaf = malloc(sizeof(LeafData)); c->leaf = malloc(sizeof(LeafData));
for (i = 0; i < MAXP; ++i) c->leaf->n = 0;
c->leaf->particle[i] = -1;
} }
Cell *cellNew(const double *xmin, const double *xmax) { Cell *cellNew(const double *xmin, const double *xmax) {
...@@ -75,7 +75,7 @@ static int intersect(const double *amin, const double *amax, const double *bmin, ...@@ -75,7 +75,7 @@ static int intersect(const double *amin, const double *amax, const double *bmin,
} }
static void leafSearch(const LeafData *l, const double *xmin, const double *xmax, int **result) { static void leafSearch(const LeafData *l, const double *xmin, const double *xmax, int **result) {
for (int i = 0; i < MAXP && l->particle[i] != -1; ++i) { for (int i = 0; i < l->n ; ++i) {
if (intersect(xmin, xmax, l->xmin[i], l->xmax[i])) if (intersect(xmin, xmax, l->xmin[i], l->xmax[i]))
vectorInsertUnique(result, l->particle[i]); vectorInsertUnique(result, l->particle[i]);
} }
...@@ -115,7 +115,7 @@ void cellSplit(Cell *c) { ...@@ -115,7 +115,7 @@ void cellSplit(Cell *c) {
} }
#endif #endif
for (i = 0; i < NSUB; ++i) for (i = 0; i < NSUB; ++i)
for (j = 0; j < MAXP; ++j) for (j = 0; j < c->leaf->n; ++j)
cellAdd(&c->sub[i], c->leaf->xmin[j], c->leaf->xmax[j], c->leaf->particle[j], NULL); cellAdd(&c->sub[i], c->leaf->xmin[j], c->leaf->xmax[j], c->leaf->particle[j], NULL);
free(c->leaf); free(c->leaf);
c->leaf = NULL; c->leaf = NULL;
...@@ -128,15 +128,15 @@ void cellAdd(Cell *c, const double *xmin, const double *xmax, int id, int **resu ...@@ -128,15 +128,15 @@ void cellAdd(Cell *c, const double *xmin, const double *xmax, int id, int **resu
if (c->leaf) { if (c->leaf) {
if (result) if (result)
leafSearch(c->leaf, xmin, xmax, result); leafSearch(c->leaf, xmin, xmax, result);
for (i = 0; i < MAXP; ++i) { if (c->leaf->n != MAXP) {
if (c->leaf->particle[i] == -1) { int i = c->leaf->n;
c->leaf->particle[i] = id; c->leaf->particle[i] = id;
for (j = 0; j < DIMENSION; ++j) { for (j = 0; j < DIMENSION; ++j) {
c->leaf->xmin[i][j] = xmin[j]; c->leaf->xmin[i][j] = xmin[j];
c->leaf->xmax[i][j] = xmax[j]; c->leaf->xmax[i][j] = xmax[j];
}
return;
} }
c->leaf->n++;
return;
} }
cellSplit(c); cellSplit(c);
} }
......
...@@ -10,7 +10,7 @@ ...@@ -10,7 +10,7 @@
#define n_fields (DIMENSION+1) #define n_fields (DIMENSION+1)
#define newton_max_it 10 #define newton_max_it 10
#define newton_atol 1e-6 #define newton_atol 2.5e-6
#define newton_rtol 1e-5 #define newton_rtol 1e-5
#if DIMENSION==2 #if DIMENSION==2
...@@ -467,6 +467,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, ...@@ -467,6 +467,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
Mesh *mesh = mesh_load_msh(mesh_file_name); Mesh *mesh = mesh_load_msh(mesh_file_name);
if (!mesh) if (!mesh)
return NULL; return NULL;
problem->mesh_tree = mesh_tree_create(mesh);
problem->rho = rho; problem->rho = rho;
problem->mu = mu; problem->mu = mu;
...@@ -543,6 +544,7 @@ void fluid_problem_free(FluidProblem *problem) ...@@ -543,6 +544,7 @@ void fluid_problem_free(FluidProblem *problem)
for (int i = 0; i < problem->n_strong_boundaries; ++i) for (int i = 0; i < problem->n_strong_boundaries; ++i)
free(problem->strong_boundaries[i].tag); free(problem->strong_boundaries[i].tag);
free(problem->strong_boundaries); free(problem->strong_boundaries);
mesh_tree_free(problem->mesh_tree);
mesh_free(problem->mesh); mesh_free(problem->mesh);
} }
...@@ -743,7 +745,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, ...@@ -743,7 +745,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
clock_get(&toc); clock_get(&toc);
printf("Time newx : %.6e \n",clock_diff(&tic,&toc)); printf("Time newx : %.6e \n",clock_diff(&tic,&toc));
clock_get(&tic); clock_get(&tic);
mesh_particles_to_mesh(problem->mesh, new_mesh->n_nodes, newx, new_eid, new_xi); mesh_tree_particle_to_mesh(problem->mesh_tree,new_mesh->n_nodes, newx, new_eid, new_xi);
clock_get(&toc); clock_get(&toc);
printf("Time find new point in old mesh : %.6e \n",clock_diff(&tic,&toc)); printf("Time find new point in old mesh : %.6e \n",clock_diff(&tic,&toc));
for (int i = 0; i < new_mesh->n_nodes; ++i) { for (int i = 0; i < new_mesh->n_nodes; ++i) {
...@@ -778,9 +780,11 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, ...@@ -778,9 +780,11 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
} }
mesh_free(problem->mesh); mesh_free(problem->mesh);
problem->mesh = new_mesh; problem->mesh = new_mesh;
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(problem->mesh);
for (int i = 0; i < problem->n_particles; ++i) for (int i = 0; i < problem->n_particles; ++i)
problem->particle_element_id[i] = -1; problem->particle_element_id[i] = -1;
mesh_particles_to_mesh(problem->mesh, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw); mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_node_volume(problem); fluid_problem_compute_node_volume(problem);
#if DIMENSION == 2 #if DIMENSION == 2
for (int i = 0; i < problem->mesh->n_nodes; ++i){ for (int i = 0; i < problem->mesh->n_nodes; ++i){
...@@ -837,7 +841,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou ...@@ -837,7 +841,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
} }
} }
clock_get(&tic); clock_get(&tic);
mesh_particles_to_mesh(problem->mesh, n, position, problem->particle_element_id, problem->particle_uvw); mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
clock_get(&toc); clock_get(&toc);
printf("Time mesh_particles_to_mesh : %.6e \n",clock_diff(&tic,&toc)); printf("Time mesh_particles_to_mesh : %.6e \n",clock_diff(&tic,&toc));
clock_get(&tic); clock_get(&tic);
......
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#define FLUID_PROBLEM_H #define FLUID_PROBLEM_H
#include "mesh.h" #include "mesh.h"
#include "mesh_find.h"
#include "hxt_linear_system.h" #include "hxt_linear_system.h"
typedef struct { typedef struct {
...@@ -17,6 +18,7 @@ typedef struct { ...@@ -17,6 +18,7 @@ typedef struct {
double alpha; double alpha;
double g; double g;
Mesh *mesh; Mesh *mesh;
MeshTree *mesh_tree;
double *porosity; double *porosity;
double *old_porosity; double *old_porosity;
double *solution; double *solution;
......
...@@ -3,6 +3,13 @@ ...@@ -3,6 +3,13 @@
#include "quadtree.h" #include "quadtree.h"
#include "mesh.h" #include "mesh.h"
#include "vector.h" #include "vector.h"
#include "mesh_find.h"
struct MeshTree_{
Cell *tree;
double *elements;
Mesh *mesh;
};
static void _bboxtol(double *bbmin, double *bbmax) { static void _bboxtol(double *bbmin, double *bbmax) {
double lmax = 0; double lmax = 0;
...@@ -16,8 +23,7 @@ static void _bboxtol(double *bbmin, double *bbmax) { ...@@ -16,8 +23,7 @@ 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]) {
double volTet(double x0[3], double x1[3], double x2[3], double x3[3]) {
// x1 * (y2 * (z3 - z4) - z2 * (y3 - y4) + y3 * z4 - y4 * z3) // x1 * (y2 * (z3 - z4) - z2 * (y3 - y4) + y3 * z4 - y4 * z3)
//- y1 * (x2 * (z3 - z4) - z2 * (x3 - x4) + x3 * z4 - x4 * z3) //- y1 * (x2 * (z3 - z4) - z2 * (x3 - x4) + x3 * z4 - x4 * z3)
//+ z1 * (x2 * (y3 - y4) - y2 * (x3 - x4) + x3 * y4 - x4 * y3) //+ z1 * (x2 * (y3 - y4) - y2 * (x3 - x4) + x3 * y4 - x4 * y3)
...@@ -27,48 +33,46 @@ double volTet(double x0[3], double x1[3], double x2[3], double x3[3]) { ...@@ -27,48 +33,46 @@ double volTet(double x0[3], double x1[3], double x2[3], double x3[3]) {
+ x0[2] * (x1[0] * (x2[1] - x3[1]) - x1[1] * (x2[0] - x3[0]) + x2[0] * x3[1] - x3[0] * x2[1]) + 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])); - (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);
if (v0 / v < -1e-8)
return 0;
double v1 = volTet(x0, x, x2, x3);
if (v1 / v < -1e-8)
return 0;
double v2 = volTet(x0, x1, x, x3);
if (v2 / v < -1e-8)
return 0;
double v3 = volTet(x0, x1, x2, x);
if (v3 / v < -1e-8)
return 0;
xi[0] = v1 / v;
xi[1] = v2 / v;
xi[2] = v3 / v;
return 1;
}
void particle_to_mesh(size_t nParticles, double *particles, int nElements, double *elements, int *elid, double *Xi) void mesh_tree_particle_to_mesh(MeshTree *mt, size_t nParticles, double *particles, int *elid, double *Xi)
{ {
double bbmin[DIMENSION], bbmax[DIMENSION];
int N = DIMENSION + 1;
for (int i = 0; i < DIMENSION; ++i) {
bbmin[i] = elements[i];
bbmax[i] = elements[i];
}
for (int i = 0; i < N * nElements; ++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]);
}
}
_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 d = 0; d < DIMENSION; ++d) {
amin[d] = el[d];
amax[d] = el[d];
}
for (int v = 1; v < N; ++v) {//only simplices
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]);
}
}
_bboxtol(amin, amax);
cellAdd(tree, amin, amax, i, NULL);
}
int *found = NULL; int *found = NULL;
vectorPushN(&found, 100); vectorPushN(&found, 100);
vectorClear(found); vectorClear(found);
int N = DIMENSION + 1;
const double *elements = mt->elements;
for (size_t i = 0; i < nParticles; ++i) { for (size_t i = 0; i < nParticles; ++i) {
double *x = &particles[i * DIMENSION]; double *x = &particles[i * DIMENSION];
if (elid[i]!=-1) { if (elid[i]!=-1) {
double *el = &elements[elid[i] * N * DIMENSION]; const double *el = &elements[elid[i] * N * DIMENSION];
#if DIMENSION == 2 #if DIMENSION == 2
double *X[3] = {el, el + DIMENSION, el + 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 = 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 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 DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
...@@ -82,33 +86,19 @@ void particle_to_mesh(size_t nParticles, double *particles, int nElements, doubl ...@@ -82,33 +86,19 @@ void particle_to_mesh(size_t nParticles, double *particles, int nElements, doubl
else else
elid[i] = -1; elid[i] = -1;
#else #else
double *X[4] = {el, el + DIMENSION, el + DIMENSION * 2, el + DIMENSION * 3}; const double *X[4] = {el, el + DIMENSION, el + DIMENSION * 2, el + DIMENSION * 3};
double v = volTet(X[0], X[1], X[2], X[3]); double xi[3];
double v0 = volTet(x, X[1], X[2], X[3]); if (!isInsideTet(x,el,el+3,el+6,el+9,Xi+i*DIMENSION))
if (v0 / v < -1e-8) elid[i] = -1;
elid[i] = -1;
double v1 = volTet(X[0], x, X[2], X[3]);
if (v1 / v < -1e-8)
elid[i] = -1 ;
double v2 = volTet(X[0], X[1], x, X[3]);
if (v2 / v < -1e-8)
elid[i] = -1;
double v3 = volTet(X[0], X[1], X[2], x);
if (v3 / v < -1e-8)
elid[i] = -1;
Xi[i * 3 + 0] = v1 / v;
Xi[i * 3 + 1] = v2 / v;
Xi[i * 3 + 2] = v3 / v;
#endif #endif
} }
if (elid[i] != -1) if (elid[i] != -1)
continue; continue;
cellSearch(tree, x, x, &found); cellSearch(mt->tree, x, x, &found);
for (size_t j = 0; j < vectorSize(found); ++j) { for (size_t j = 0; j < vectorSize(found); ++j) {
double *el = &elements[found[j] * N * DIMENSION]; const double *el = &elements[found[j] * N * DIMENSION];
if (DIMENSION == 2) { if (DIMENSION == 2) {
double *X[3] = {el, el + DIMENSION, el + 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 = 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 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 DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
...@@ -123,25 +113,10 @@ void particle_to_mesh(size_t nParticles, double *particles, int nElements, doubl ...@@ -123,25 +113,10 @@ void particle_to_mesh(size_t nParticles, double *particles, int nElements, doubl
} }
} }
else { else {
double *X[4] = {el, el + DIMENSION, el + DIMENSION * 2, el + DIMENSION * 3}; if (isInsideTet(x,el,el+3,el+6,el+9,Xi+i*DIMENSION)){
double v = volTet(X[0], X[1], X[2], X[3]); elid[i] = found[j];
double v0 = volTet(x, X[1], X[2], X[3]); break;
if (v0 / v < -1e-8) }
continue;
double v1 = volTet(X[0], x, X[2], X[3]);
if (v1 / v < -1e-8)
continue;
double v2 = volTet(X[0], X[1], x, X[3]);
if (v2 / v < -1e-8)
continue;
double v3 = volTet(X[0], X[1], X[2], x);
if (v3 / v < -1e-8)
continue;
Xi[i * 3 + 0] = v1 / v;
Xi[i * 3 + 1] = v2 / v;
Xi[i * 3 + 2] = v3 / v;
elid[i] = found[j];
break;
} }
} }
/*if (elid[i] == -1) { /*if (elid[i] == -1) {
...@@ -169,23 +144,24 @@ void particle_to_mesh(size_t nParticles, double *particles, int nElements, doubl ...@@ -169,23 +144,24 @@ void particle_to_mesh(size_t nParticles, double *particles, int nElements, doubl
//exit(1); //exit(1);
} }
//} //}
vectorClear(found);
} }
vectorFree(found); vectorFree(found);
cellFree(tree);
} }
void mesh_particles_to_mesh(const Mesh *mesh, size_t nParticles, double *particles, int *elid, double *xi)
{ MeshTree *mesh_tree_create(Mesh *mesh) {
#if DIMENSION == 2 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); double *elements = malloc(sizeof(double)*6*mesh->n_elements);
for (int i = 0; i < mesh->n_elements*3; ++i) { for (int i = 0; i < nElements*3; ++i) {
int iv = mesh->elements[i]; int iv = mesh->elements[i];
elements[i*2+0] = mesh->x[iv*3+0]; elements[i*2+0] = mesh->x[iv*3+0];
elements[i*2+1] = mesh->x[iv*3+1]; elements[i*2+1] = mesh->x[iv*3+1];
} }
particle_to_mesh(nParticles,particles,mesh->n_elements,elements,elid,xi); #else
free(elements);
#else
double *elements = malloc(sizeof(double)*12*mesh->n_elements); double *elements = malloc(sizeof(double)*12*mesh->n_elements);
for (int i = 0; i < mesh->n_elements*4; ++i) { for (int i = 0; i < mesh->n_elements*4; ++i) {
int iv = mesh->elements[i]; int iv = mesh->elements[i];
...@@ -193,7 +169,45 @@ void mesh_particles_to_mesh(const Mesh *mesh, size_t nParticles, double *particl ...@@ -193,7 +169,45 @@ void mesh_particles_to_mesh(const Mesh *mesh, size_t nParticles, double *particl
elements[i*3+1] = mesh->x[iv*3+1]; elements[i*3+1] = mesh->x[iv*3+1];
elements[i*3+2] = mesh->x[iv*3+2]; elements[i*3+2] = mesh->x[iv*3+2];
} }
particle_to_mesh(nParticles,particles,mesh->n_elements,elements,elid,xi); #endif
free(elements); mt->elements = elements;
#endif double bbmin[DIMENSION], bbmax[DIMENSION];
int N = DIMENSION + 1;
for (int i = 0; i < DIMENSION; ++i) {
bbmin[i] = elements[i];
bbmax[i] = elements[i];
}
for (int i = 0; i < N * nElements; ++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]);
}
}
_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 d = 0; d < DIMENSION; ++d) {
amin[d] = el[d];
amax[d] = el[d];
}
for (int v = 1; v < N; ++v) {//only simplices
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]);
}
}
_bboxtol(amin, amax);
cellAdd(tree, amin, amax, i, NULL);
}
mt->tree = tree;
return mt;
}
void mesh_tree_free(MeshTree *mt) {
cellFree(mt->tree);
free(mt->elements);
free(mt);
} }
#ifndef MESH_FIND_H #ifndef MESH_FIND_H
#define MESH_FIND_H #define MESH_FIND_H
#include "mesh.h" #include "mesh.h"
void particles_to_mesh(size_t nParticles, double *particles, int nElements, double *elements, int *elid, double *Xi); typedef struct MeshTree_ MeshTree;
void mesh_particles_to_mesh(const Mesh *mesh, size_t n_particles, double *particles, int *elid, double *xi); MeshTree *mesh_tree_create(Mesh *mesh);
void mesh_tree_particle_to_mesh(MeshTree *mt, size_t nParticles, double *particles, int *elid, double *Xi);
void mesh_tree_free(MeshTree *mt);
#endif #endif
Supports Markdown
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