[Bf-blender-cvs] [2b2b3ea16b1] soc-2022-soft-bodies: Fixed a few bugs. Added floating point error handling. Added view for multiple faces per tet for better visualization. Added support for removing external tet removal using BVH tree
Aarnav Dhanuka
noreply at git.blender.org
Mon Jul 11 10:19:46 CEST 2022
Commit: 2b2b3ea16b119cdf37f726d030d82c5c6fa1cec6
Author: Aarnav Dhanuka
Date: Mon Jul 11 13:49:44 2022 +0530
Branches: soc-2022-soft-bodies
https://developer.blender.org/rB2b2b3ea16b119cdf37f726d030d82c5c6fa1cec6
Fixed a few bugs. Added floating point error handling. Added view for multiple faces per tet for better visualization. Added support for removing external tet removal using BVH tree
===================================================================
M source/blender/modifiers/intern/MOD_weld.cc
===================================================================
diff --git a/source/blender/modifiers/intern/MOD_weld.cc b/source/blender/modifiers/intern/MOD_weld.cc
index ce1b02244bb..b6e70dda538 100644
--- a/source/blender/modifiers/intern/MOD_weld.cc
+++ b/source/blender/modifiers/intern/MOD_weld.cc
@@ -38,10 +38,7 @@
#include "DNA_modifier_types.h"
#include "DNA_screen_types.h"
-#ifdef USE_BVHTREEKDOP
-# include "BKE_bvhutils.h"
-#endif
-
+#include "BKE_bvhutils.h"
#include "BKE_context.h"
#include "BKE_deform.h"
#include "BKE_modifier.h"
@@ -71,15 +68,40 @@ using namespace blender::math;
Vector<Vector<int>> tetFaces({{2,1,0}, {0,1,3}, {1,2,3}, {2,0,3}});
-float inf = 1000000000000;// C has no FLOAT_MAX :/
+float directions[6][3] = {{1,0,0}, {-1,0,0}, {0,1,0}, {0,-1,0}, {0,0,1}, {0,0,-1}};
+
+float inf = FLT_MAX;// C has no FLOAT_MAX :/
+float eps = 1e-4;
static float randomEps(){
- // Why are we using a random epsilon value?
- return 0;
+ float eps = 1e-6 ;
+ return -eps + 2.0 * (static_cast <float> (rand()) / static_cast <float> (RAND_MAX)) * eps;
}
-static bool isInside(Vector<float> vert){
- return true;
+static bool isInside(float3 vert, BVHTreeFromMesh *treedata){
+ int count = 0;
+ float min_dist = 0.0f;
+ for(auto dir : directions){
+ float radius = 0.0f;
+
+ float max_length = FLT_MAX;
+
+ BVHTreeRayHit rayhit = {0};
+ rayhit.index = -1;
+ rayhit.dist = max_length;
+ BLI_bvhtree_ray_cast(treedata->tree, vert, dir, radius, &rayhit, treedata->raycast_callback, treedata);
+
+ if (rayhit.index != -1 && rayhit.dist <= max_length) {
+ if(dot_v3v3(rayhit.no, dir) > eps){
+ count++;
+ }
+
+ if((min_dist > eps) && (min_dist - rayhit.dist) > eps)
+ return false;
+ }
+ }
+
+ return count > 3;
}
// Used to populate a tets face normals and planesD
@@ -101,6 +123,31 @@ static void setTetProperties(Vector<float3> &verts,
}
}
+static float3 getCircumCenter(float3 p0, float3 p1, float3 p2, float3 p3){
+ float3 b = p1 - p0;
+ float3 c = p2 - p0;
+ float3 d = p3 - p0;
+
+ float det = 2.0 * (b.x*(c.y*d.z - c.z*d.y) - b.y*(c.x*d.z - c.z*d.x) + b.z*(c.x*d.y - c.y*d.x));
+ if (det <= eps){
+ return p0;
+ }
+ else{
+ float3 v = cross(c, d)*dot(b, b) + cross(d, b)*dot(c, c) + cross(b, c)*dot(d, d);
+ v /= det;
+ return p0 + v;
+ }
+
+}
+
+// bool isSameSide(float3 p0, float3 p1, float3 p2, float3 p4, float3 vert){
+
+// }
+
+// bool isInsideTet(float3 p0, float3 p1, float3 p2, float3 p3, float3 vert){
+
+// }
+
static int findContainingTet(Vector<float3> &verts,
Vector<int> &tetVertId,
Vector<int> &tetFaceNeighbors,
@@ -110,85 +157,121 @@ static int findContainingTet(Vector<float3> &verts,
Vector<int> &tetMarks,
float3 currVert){
- bool found = false;
- int tetNr = 0;
+ /* --------------------
+ Matthias Method
+ ----------------------*/
- while(tetNr < tetVertId.size()/4 && tetVertId[4*tetNr]<0)
- tetNr++;
+ // bool found = false;
+ // int tetNr = 0;
- float3 center(0.0, 0.0, 0.0);
- while(!found){
- if(tetNr < 0 || tetMarks[tetNr] == tetMarkId){
- break;
- }
-
- tetMarks[tetNr] = tetMarkId;
+ // while(tetNr < tetVertId.size()/4 && tetVertId[4*tetNr]<0)
+ // tetNr++;
- center = {0.0,0.0,0.0};
- for(int i = 0; i<4; i++){
- center += verts[tetVertId[4*tetNr + i]];
- }
- center*=0.25;
+ // float3 center(0.0, 0.0, 0.0);
+ // while(!found){
+ // if(tetNr < 0 || tetMarks[tetNr] == tetMarkId){
+ // break;
+ // }
+
+ // tetMarks[tetNr] = tetMarkId;
+
+ // center = {0.0,0.0,0.0};
+ // for(int i = 0; i<4; i++){
+ // center += verts[tetVertId[4*tetNr + i]];
+ // }
+ // center*=0.25;
+
+ // float minT = inf; //
+ // int minFaceNr = -1;
+
+ // for(int i = 0; i<4; i++){
+ // float3 normal = faceNormals[4*tetNr + i];
+ // float d = planesD[4*tetNr + i];
+
+ // float hp = dot(normal, currVert) - d;
+ // float hc = dot(normal, center) - d;
+
+ // float t = hp - hc;
+ // if(t <= eps)
+ // continue;
+
+ // t = -hc/t;
+
+ // if((t >= eps) && ((t - minT) > eps)){
+ // minT = t;
+ // minFaceNr = i;
+ // }
+ // }
+
+ // if(minT >= 1.0){
+ // found = true;
+ // }
+ // else{
+ // tetNr = tetFaceNeighbors[4*tetNr + minFaceNr];
+ // }
+ // }
- float minT = inf; //
- int minFaceNr = -1;
+ // if(found)
+ // return tetNr;
+ // else
+ // return -1;
- for(int i = 0; i<4; i++){
- float3 normal = faceNormals[4*tetNr + i];
- float d = planesD[4*tetNr + i];
+ /* --------------------
+ Brute force finding first violating tet
+ ----------------------*/
- float hp = dot(normal, currVert) - d;
- float hc = dot(normal, center) - d;
+ for(int currTet = 0; currTet<tetVertId.size()/4; currTet++){
+ if(tetVertId[4*currTet] == -1)
+ continue;
+
+ float3 p0 = verts[tetVertId[4*currTet + 0]];
+ float3 p1 = verts[tetVertId[4*currTet + 1]];
+ float3 p2 = verts[tetVertId[4*currTet + 2]];
+ float3 p3 = verts[tetVertId[4*currTet + 3]];
- float t = hp - hc;
- if(t == 0.0)
- continue;
+ float3 circumCenter = getCircumCenter(p0, p1, p2, p3);
+ float circumRadius = length(p0 - circumCenter);
+
+ if((circumRadius - length(currVert - circumCenter)) >= eps){
+ return currTet;
+ }
+ }
- t = -hc/t;
+ /* -----------------------
+ My method, checks if point is inside tet, if not finds the face that is connecting to the point and center
+ ------------------------*/
- if(t >= 0.0 && t < minT){
- minT = t;
- minFaceNr = i;
- }
- }
+ // bool found = false;
+ // int tetNr = 0;
+ // while(tetNr < tetVertId.size()/4 && tetVertId[4*tetNr] < 0)
+ // tetNr++;
- if(minT >= 1.0){
- found = true;
- }
- else{
- tetNr = tetFaceNeighbors[4*tetNr + minFaceNr];
- }
- }
+ // while(!found){
+ // float3 p0 = verts[tetVertId[4*currTet + 0]];
+ // float3 p1 = verts[tetVertId[4*currTet + 1]];
+ // float3 p2 = verts[tetVertId[4*currTet + 2]];
+ // float3 p3 = verts[tetVertId[4*currTet + 3]];
+
+ // if(isInsideTet(p0, p1, p2, p3, currVert))
+ // return tetNr;
+
+ // for(int i = 0; i<4; i++)
+ // }
- if(found)
- return tetNr;
- else
- return -1;
+ return -1;
}
-static float3 getCircumCenter(float3 p0, float3 p1, float3 p2, float3 p3){
- float3 b = p1 - p0;
- float3 c = p2 - p0;
- float3 d = p3 - p0;
+/*
+The basic assumption employed is that 2 violating tets cannot not be neighbors.
+Simple BFS approach that checks neighbors of all violating tets and adds them to stack
- float det = 2.0 * (b.x*(c.y*d.z - c.z*d.y) - b.y*(c.x*d.z - c.z*d.x) + b.z*(c.x*d.y - c.y*d.x));
- if (det == 0.0){
- return p0;
- }
- else{
- float3 v = cross(c, d)*dot(b, b) + cross(d, b)*dot(c, c) + cross(b, c)*dot(d, d);
- v /= det;
- return p0 + v;
- }
-
-}
+If a violating tet shares a face with a non-violating tet, that's a boundary face and for each violating tet, list of its boundary faces is returned
-// Remove unused params
+Note - BFS can be written better
+*/
static Vector<std::pair<int, Vector<int>>> getViolatingTets(Vector<float3> &verts,
Vector<int> &tetVertId,
Vector<int> &tetFaceNeighbors,
- Vector<float3> &faceNormals,
- Vector<float> &planesD,
int tetMarkId,
Vector<int> &tetMarks,
float3 currVert,
@@ -198,13 +281,15 @@ static Vector<std::pair<int, Vector<int>>> getViolatingTets(Vector<float3> &vert
Vector< std::pair<int,Vector<int>> > violatingTets;
Vector<int> stack;
- stack.append(containingTetNr);
- tetMarks[containingTetNr] = tetMarkId;
+ stack.append(containingTetNr);
while(stack.size()){
int currTet = stack.last();
stack.remove_last();
+ if(tetMarks[currTet] == tetMarkId){
+ continue;
+ }
tetMarks[currTet] = tetMarkId;
Vector<int> currTetBorderFaces;
@@ -219,17 +304,17 @@ static Vector<std::pair<int, Vector<int>>> getViolatingTets(Vector<float3> &vert
continue;
}
- float3 p0 = verts[tetVertId[4*currTet + 0]];
- float3 p1 = verts[tetVertId[4*currTet + 1]];
- float3 p2 = verts[tetVertId[4*currTet + 2]];
- float3 p3 = verts[tetVertId[4*currTet + 3]];
+ float3 p0 = verts[tetVertId[4*neighborTet + 0]];
+ float3 p1 = verts[tetVertId[4*neighborTet + 1]];
+ float3 p2 = verts[tetVertId[4*neighborTet + 2]];
+ float3 p3 = verts[tetVertId[4*neighborTet + 3]];
float3 circumCenter = getCircumCenter(p0, p1, p2, p3);
float circumRadius = length(p0 - circumCenter);
- if(length(currVert - circumCenter) <= circumRadius){
+ if((circumRadius - length(currVert - circumCenter)) >= eps){
stack.append(neighborTet);
- tetMarks[neighborTet] = tetMarkId;
+ // tetMarks[neighborTet] = tetMarkId;
}
else{
currTetBorderFaces.append(i);
@@ -238,21 +323,10 @@ static Vector<std::pair<int, Vector<int>>> getViolatingTets(Vector<float3> &vert
violatingTets.append({currTet, currTetBorderFaces});
}
- // for(int currTet = 0; currTet<tetVertId.size()/4; currTet++){
- // float3 p0 = verts[tetVertId[4*currTet + 0]];
- // float3 p1 = verts[tetVertId[4*currTet + 1]];
- // float3 p2 = verts[tetVertId[4*currTet + 2]];
- // float3 p3 = verts[tetVertId[4*currTet + 3]];
-
- // float3 circumCenter = getCircumCenter(p0, p1, p2, p3);
- // float circumRadius = length(p0 - circumCenter);
- // if()
- // }
-
return violatingTets;
}
-static Vector<int> createTets(Vector<float3> verts, /*tree, */ float minTetQuality){
+static Vector<int> createTets(Vector<float3> verts, BVHTreeFromMesh *treedata, float minTetQuality){
Vector<int> tetVertId; // Stores indices of vertex that form a tet. Every tet is stored as 4 indices
Vector<int> tetFaceNeighbors; // Stores index of tet that shares
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list