Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add new refine surface, bounding box, and point refine methods. #297

Merged
merged 9 commits into from
Dec 2, 2022
Merged

Add new refine surface, bounding box, and point refine methods. #297

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
71eb2f8
Add new refine surface, bounding box, and point refine methods.
jcapriot Nov 23, 2022
d2a4fd9
remove refine_tree_xyz in the examples.
jcapriot Nov 23, 2022
8f8a8e5
fix deprecation warning
jcapriot Nov 23, 2022
64a699c
fix syntax
jcapriot Nov 23, 2022
c7e5bfa
fix working error test
jcapriot Nov 23, 2022
7f3a956
add tests for new functions
jcapriot Nov 23, 2022
597a37f
modify test to hit last line of code
jcapriot Nov 23, 2022
4414192
refine updates
jcapriot Dec 1, 2022
b3d059c
documentation updates
jcapriot Dec 1, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
177 changes: 177 additions & 0 deletions discretize/_extensions/tree.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -685,6 +685,163 @@ void Cell::refine_triangle(
}
}

void Cell::refine_vert_triang_prism(
node_map_t& nodes,
double* x0, double* x1, double* x2, double h,
double* e0, double* e1, double* e2, double* t_norm,
int_t p_level, double *xs, double *ys, double* zs, bool diag_balance
){
// Return If I'm at max_level or p_level
if (level >= p_level || level == max_level){
return;
}
// check all the AABB faces
double v0[3], v1[3], v2[3], half[3];
double vmin, vmax;
double p0, p1, p2, p3, pmin, pmax, rad;
for(int_t i=0; i < n_dim; ++i){
v0[i] = x0[i] - location[i];
v1[i] = x1[i] - location[i];
vmin = std::min(v0[i], v1[i]);
vmax = std::max(v0[i], v1[i]);
v2[i] = x2[i] - location[i];
vmin = std::min(vmin, v2[i]);
vmax = std::max(vmax, v2[i]);
if(i == 2){
vmax += h;
}
half[i] = location[i] - points[0]->location[i];

// Bounding box check
if (vmin > half[i] || vmax < -half[i]){
return;
}
}
// first do the 3 edge cross tests that apply in 2D and 3D

// edge 0 cross z_hat
//p0 = e0[1] * v0[0] - e0[0] * v0[1];
p1 = e0[1] * v1[0] - e0[0] * v1[1];
p2 = e0[1] * v2[0] - e0[0] * v2[1];
pmin = std::min(p1, p2);
pmax = std::max(p1, p2);
rad = std::abs(e0[1]) * half[0] + std::abs(e0[0]) * half[1];
if (pmin > rad || pmax < -rad){
return;
}

// edge 1 cross z_hat
p0 = e1[1] * v0[0] - e1[0] * v0[1];
p1 = e1[1] * v1[0] - e1[0] * v1[1];
//p2 = e1[1] * v2[0] - e1[0] * v2[1];
pmin = std::min(p0, p1);
pmax = std::max(p0, p1);
rad = std::abs(e1[1]) * half[0] + std::abs(e1[0]) * half[1];
if (pmin > rad || pmax < -rad){
return;
}

// edge 2 cross z_hat
//p0 = e2[1] * v0[0] - e2[0] * v0[1];
p1 = e2[1] * v1[0] - e2[0] * v1[1];
p2 = e2[1] * v2[0] - e2[0] * v2[1];
pmin = std::min(p1, p2);
pmax = std::max(p1, p2);
rad = std::abs(e2[1]) * half[0] + std::abs(e2[0]) * half[1];
if (pmin > rad || pmax < -rad){
return;
}

// edge 0 cross x_hat
p0 = e0[2] * v0[1] - e0[1] * v0[2];
p1 = e0[2] * v0[1] - e0[1] * (v0[2] + h);
p2 = e0[2] * v2[1] - e0[1] * v2[2];
p3 = e0[2] * v2[1] - e0[1] * (v2[2] + h);
pmin = std::min(std::min(std::min(p0, p1), p2), p3);
pmax = std::max(std::max(std::max(p0, p1), p2), p3);
rad = std::abs(e0[2]) * half[1] + std::abs(e0[1]) * half[2];
if (pmin > rad || pmax < -rad){
return;
}
// edge 0 cross y_hat
p0 = -e0[2] * v0[0] + e0[0] * v0[2];
p1 = -e0[2] * v0[0] + e0[0] * (v0[2] + h);
p2 = -e0[2] * v2[0] + e0[0] * v2[2];
p3 = -e0[2] * v2[0] + e0[0] * (v2[2] + h);
pmin = std::min(std::min(std::min(p0, p1), p2), p3);
pmax = std::max(std::max(std::max(p0, p1), p2), p3);
rad = std::abs(e0[2]) * half[0] + std::abs(e0[0]) * half[2];
if (pmin > rad || pmax < -rad){
return;
}
// edge 1 cross x_hat
p0 = e1[2] * v0[1] - e1[1] * v0[2];
p1 = e1[2] * v0[1] - e1[1] * (v0[2] + h);
p2 = e1[2] * v2[1] - e1[1] * v2[2];
p3 = e1[2] * v2[1] - e1[1] * (v2[2] + h);
pmin = std::min(std::min(std::min(p0, p1), p2), p3);
pmax = std::max(std::max(std::max(p0, p1), p2), p3);
rad = std::abs(e1[2]) * half[1] + std::abs(e1[1]) * half[2];
if (pmin > rad || pmax < -rad){
return;
}
// edge 1 cross y_hat
p0 = -e1[2] * v0[0] + e1[0] * v0[2];
p1 = -e1[2] * v0[0] + e1[0] * (v0[2] + h);
p2 = -e1[2] * v2[0] + e1[0] * v2[2];
p3 = -e1[2] * v2[0] + e1[0] * (v2[2] + h);
pmin = std::min(std::min(std::min(p0, p1), p2), p3);
pmax = std::max(std::max(std::max(p0, p1), p2), p3);
rad = std::abs(e1[2]) * half[0] + std::abs(e1[0]) * half[2];
if (pmin > rad || pmax < -rad){
return;
}
// edge 2 cross x_hat
p0 = e2[2] * v0[1] - e2[1] * v0[2];
p1 = e2[2] * v0[1] - e2[1] * (v0[2] + h);
p2 = e2[2] * v1[1] - e2[1] * v1[2];
p3 = e2[2] * v1[1] - e2[1] * (v1[2] + h);
pmin = std::min(std::min(std::min(p0, p1), p2), p3);
pmax = std::max(std::max(std::max(p0, p1), p2), p3);
rad = std::abs(e2[2]) * half[1] + std::abs(e2[1]) * half[2];
if (pmin > rad || pmax < -rad){
return;
}
// edge 2 cross y_hat
p0 = -e2[2] * v0[0] + e2[0] * v0[2];
p1 = -e2[2] * v0[0] + e2[0] * (v0[2] + h);
p2 = -e2[2] * v1[0] + e2[0] * v1[2];
p3 = -e2[2] * v1[0] + e2[0] * (v1[2] + h);
pmin = std::min(std::min(std::min(p0, p1), p2), p3);
pmax = std::max(std::max(std::max(p0, p1), p2), p3);
rad = std::abs(e2[2]) * half[0] + std::abs(e2[0]) * half[2];
if (pmin > rad || pmax < -rad){
return;
}

// triangle normal axis
p0 = t_norm[0] * v0[0] + t_norm[1] * v0[1] + t_norm[2] * v0[2];
p1 = t_norm[0] * v0[0] + t_norm[1] * v0[1] + t_norm[2] * (v0[2] + h);
pmin = std::min(p0, p1);
pmax = std::max(p0, p1);
rad = std::abs(t_norm[0]) * half[0] + std::abs(t_norm[1]) * half[1] + std::abs(t_norm[2]) * half[2];
if (pmin > rad || pmax < -rad){
return;
}
// the axes defined by the three vertical prism faces
// should already be tested by the e0, e1, e2 cross z_hat tests

// If here, then I intersect the triangle!
if(is_leaf()){
divide(nodes, xs, ys, zs, true, diag_balance);
}
for(int_t i = 0; i < (1<<n_dim); ++i){
children[i]->refine_vert_triang_prism(
nodes, x0, x1, x2, h, e0, e1, e2, t_norm, p_level, xs, ys, zs, diag_balance
);
}
}

void Cell::refine_tetra(
node_map_t& nodes,
double* x0, double* x1, double* x2, double* x3,
Expand Down Expand Up @@ -1347,6 +1504,26 @@ void Tree::refine_triangle(double* x0, double* x1, double* x2, int_t p_level, bo
);
};

void Tree::refine_vert_triang_prism(double* x0, double* x1, double* x2, double h, int_t p_level, bool diagonal_balance){
double e0[3], e1[3], e2[3], t_norm[3];
for(int_t i=0; i<n_dim; ++i){
e0[i] = x1[i] - x0[i];
e1[i] = x2[i] - x1[i];
e2[i] = x2[i] - x0[i];
}
if(n_dim > 2){
t_norm[0] = e0[1] * e1[2] - e0[2] * e1[1];
t_norm[1] = e0[2] * e1[0] - e0[0] * e1[2];
t_norm[2] = e0[0] * e1[1] - e0[1] * e1[0];
}
for(int_t iz=0; iz<nz_roots; ++iz)
for(int_t iy=0; iy<ny_roots; ++iy)
for(int_t ix=0; ix<nx_roots; ++ix)
roots[iz][iy][ix]->refine_vert_triang_prism(
nodes, x0, x1, x2, h, e0, e1, e2, t_norm, p_level, xs, ys, zs, diagonal_balance
);
};

void Tree::refine_tetra(double* x0, double* x1, double* x2, double* x3, int_t p_level, bool diagonal_balance){
double t_edges[6][3];
double face_normals[4][3];
Expand Down
8 changes: 8 additions & 0 deletions discretize/_extensions/tree.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,11 @@ class Cell{
double* x0, double* x1, double* x2, double* e0, double* e1, double* e2, double* t_norm,
int_t p_level, double *xs, double *ys, double* zs, bool diag_balance=false
);
void refine_vert_triang_prism(node_map_t& nodes,
double* x0, double* x1, double* x2, double h,
double* e0, double* e1, double* e2, double* t_norm,
int_t p_level, double *xs, double *ys, double* zs, bool diag_balance=false
);
void refine_tetra(
node_map_t& nodes,
double* x0, double* x1, double* x2, double* x3,
Expand Down Expand Up @@ -179,6 +184,9 @@ class Tree{
void refine_tetra(
double* x0, double* x1, double* x2, double* x3, int_t p_level, bool diag_balance=false
);
void refine_vert_triang_prism(
double* x0, double* x1, double* x2, double h, int_t p_level, bool diagonal_balance=false
);
void number();
void finalize_lists();

Expand Down
1 change: 1 addition & 0 deletions discretize/_extensions/tree.pxd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ cdef extern from "tree.h":
void refine_box(double*, double*, int_t, bool)
void refine_line(double*, double*, int_t, bool)
void refine_triangle(double*, double*, double*, int_t, bool)
void refine_vert_triang_prism(double*, double*, double*, double, int_t, bool)
void refine_tetra(double*, double*, double*, double*, int_t, bool)
void number()
void initialize_roots()
Expand Down
Loading