Skip to content
Snippets Groups Projects
Commit e02840df authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

cleanup

parent 28d5ade2
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/GRbf.cpp
+ 78
− 234
View file @ e02840df
...@@ -86,8 +86,9 @@ static void exportParametrizedMesh(fullMatrix<double> &UV, int nbNodes) ...@@ -86,8 +86,9 @@ static void exportParametrizedMesh(fullMatrix<double> &UV, int nbNodes)
GRbf::GRbf(double sizeBox, int variableEps, int rbfFun, GRbf::GRbf(double sizeBox, int variableEps, int rbfFun,
std::map<MVertex*, SVector3> _normals, std::map<MVertex*, SVector3> _normals,
std::set<MVertex *> allNodes, std::vector<MVertex*> bcNodes, bool _isLocal) std::set<MVertex *> allNodes, std::vector<MVertex*> bcNodes, bool _isLocal)
: sBox(sizeBox), variableShapeParam(variableEps), radialFunctionIndex (rbfFun), : isLocal(_isLocal), _inUV(0), sBox(sizeBox), variableShapeParam(variableEps),
_inUV(0), isLocal(_isLocal) radialFunctionIndex (rbfFun)
{ {
#if defined (HAVE_ANN) #if defined (HAVE_ANN)
XYZkdtree = 0; XYZkdtree = 0;
...@@ -125,7 +126,6 @@ GRbf::GRbf(double sizeBox, int variableEps, int rbfFun, ...@@ -125,7 +126,6 @@ GRbf::GRbf(double sizeBox, int variableEps, int rbfFun,
normals.resize(nbNodes,3); normals.resize(nbNodes,3);
int index = 0; int index = 0;
double dist_min = 1.e6; double dist_min = 1.e6;
double dist_max = 1.e-6;
for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end(); ++itv){ for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end(); ++itv){
MVertex *v1 = *itv; MVertex *v1 = *itv;
centers(index,0) = v1->x()/sBox; centers(index,0) = v1->x()/sBox;
...@@ -223,7 +223,6 @@ void GRbf::buildOctree(double radius) ...@@ -223,7 +223,6 @@ void GRbf::buildOctree(double radius)
if (l.size() == 1) printf("*** WARNING: Found only %d sphere ! \n", (int)l.size()); if (l.size() == 1) printf("*** WARNING: Found only %d sphere ! \n", (int)l.size());
for (std::vector<void*>::iterator it = l.begin(); it != l.end(); it++) { for (std::vector<void*>::iterator it = l.begin(); it != l.end(); it++) {
Sphere *sph = (Sphere*) *it; Sphere *sph = (Sphere*) *it;
std::vector<int> &pts = nodesInSphere[i];
if (sph->index != i) nodesInSphere[i].push_back(sph->index); if (sph->index != i) nodesInSphere[i].push_back(sph->index);
} }
//printf("size node i =%d = %d \n", i , nodesInSphere[i].size()); //printf("size node i =%d = %d \n", i , nodesInSphere[i].size());
...@@ -285,7 +284,7 @@ void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs, ...@@ -285,7 +284,7 @@ void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
fullMatrix<double> nodes_in_sph(pts.size(),3); fullMatrix<double> nodes_in_sph(pts.size(),3);
fullMatrix<double> LocalOper; fullMatrix<double> LocalOper;
for (int k=0; k< pts.size(); k++){ for (unsigned int k = 0; k < pts.size(); k++){
nodes_in_sph(k,0) = cntrs(pts[k],0); nodes_in_sph(k,0) = cntrs(pts[k],0);
nodes_in_sph(k,1) = cntrs(pts[k],1); nodes_in_sph(k,1) = cntrs(pts[k],1);
nodes_in_sph(k,2) = cntrs(pts[k],2); nodes_in_sph(k,2) = cntrs(pts[k],2);
...@@ -336,6 +335,7 @@ double GRbf::evalRadialFnDer (int p, double dx, double dy, double dz, double ep) ...@@ -336,6 +335,7 @@ double GRbf::evalRadialFnDer (int p, double dx, double dy, double dz, double ep)
case 222: return ep*ep*(3+ep*ep*2*r2)/sqrt((1+ep*ep*r2)*(1+ep*ep*r2)*(1+ep*ep*r2)); case 222: return ep*ep*(3+ep*ep*2*r2)/sqrt((1+ep*ep*r2)*(1+ep*ep*r2)*(1+ep*ep*r2));
} }
} }
return 0.;
} }
fullMatrix<double> GRbf::generateRbfMat(int p, fullMatrix<double> GRbf::generateRbfMat(int p,
...@@ -477,7 +477,7 @@ void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs, ...@@ -477,7 +477,7 @@ void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs,
LocalOper.setAll(0.0); LocalOper.setAll(0.0);
for (int k=0; k< pts.size(); k++){ for (unsigned int k = 0; k < pts.size(); k++){
nodes_in_sph(k, 0) = cntrs(pts[k], 0); nodes_in_sph(k, 0) = cntrs(pts[k], 0);
nodes_in_sph(k, 1) = cntrs(pts[k], 1); nodes_in_sph(k, 1) = cntrs(pts[k], 1);
nodes_in_sph(k, 2) = cntrs(pts[k], 2); nodes_in_sph(k, 2) = cntrs(pts[k], 2);
...@@ -488,7 +488,7 @@ void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs, ...@@ -488,7 +488,7 @@ void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs,
RbfLapSurface_global_projection(nodes_in_sph,local_normals, LocalOper); RbfLapSurface_global_projection(nodes_in_sph,local_normals, LocalOper);
for (int j=0;j<pts.size() ; j++) for (unsigned int j = 0; j < pts.size(); j++)
Oper(i, pts[j]) = LocalOper(0, j); Oper(i, pts[j]) = LocalOper(0, j);
} }
} }
...@@ -614,13 +614,13 @@ void GRbf::RbfLapSurface_local_CPM_sparse(std::vector<MVertex*> &bndVertices, bo ...@@ -614,13 +614,13 @@ void GRbf::RbfLapSurface_local_CPM_sparse(std::vector<MVertex*> &bndVertices, bo
std::vector<int> &pts = nodesInSphere[i]; std::vector<int> &pts = nodesInSphere[i];
if (bndIndices.count(i) > 0) { if (bndIndices.count(i) > 0) {
sys.insertInSparsityPattern(i, i); sys.insertInSparsityPattern(i, i);
for (int j = 0; j < pts.size(); ++j) { for (unsigned int j = 0; j < pts.size(); ++j) {
sys.insertInSparsityPattern(i + numNodes, pts[j]); sys.insertInSparsityPattern(i + numNodes, pts[j]);
sys.insertInSparsityPattern(i + 2 * numNodes, pts[j]); sys.insertInSparsityPattern(i + 2 * numNodes, pts[j]);
} }
} }
else { else {
for (int j = 0; j < pts.size(); ++j) { for (unsigned int j = 0; j < pts.size(); ++j) {
sys.insertInSparsityPattern(i, pts[j]); sys.insertInSparsityPattern(i, pts[j]);
sys.insertInSparsityPattern(i + numNodes, pts[j]); sys.insertInSparsityPattern(i + numNodes, pts[j]);
sys.insertInSparsityPattern(i + 2 * numNodes, pts[j]); sys.insertInSparsityPattern(i + 2 * numNodes, pts[j]);
...@@ -833,8 +833,8 @@ void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs, ...@@ -833,8 +833,8 @@ void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs,
for (int j = 0; j < nnTot; ++j){ for (int j = 0; j < nnTot; ++j){
Oper(i,j) = PLap(i,j); Oper(i,j) = PLap(i,j);
double del = (i == j) ? -1.0: 0.0; double del = (i == j) ? -1.0: 0.0;
double pos1 = (i+numNodes == j) ? 1.0: 0.0; //double pos1 = (i+numNodes == j) ? 1.0: 0.0;
double pos2 = (i+2*numNodes == j) ? 1.0: 0.0; //double pos2 = (i+2*numNodes == j) ? 1.0: 0.0;
Oper(i+numNodes,j) = del +Ix2extX(i,j);//+ pos1; // Oper(i+numNodes,j) = del +Ix2extX(i,j);//+ pos1; //
Oper(i+2*numNodes,j) = del + Ix2extX(i+numNodes,j); //+ pos2; // Oper(i+2*numNodes,j) = del + Ix2extX(i+numNodes,j); //+ pos2; //
} }
...@@ -854,7 +854,7 @@ void GRbf::solveHarmonicMap_sparse(linearSystem<double> &sys, int numNodes, ...@@ -854,7 +854,7 @@ void GRbf::solveHarmonicMap_sparse(linearSystem<double> &sys, int numNodes,
for (int j = 0; j < 2; ++j) { for (int j = 0; j < 2; ++j) {
sys.zeroRightHandSide(); sys.zeroRightHandSide();
//UNIT CIRCLE //UNIT CIRCLE
for (int i = 0; i < bcNodes.size(); i++) { for (unsigned int i = 0; i < bcNodes.size(); i++) {
std::set<MVertex *>::iterator itN = myNodes.find(bcNodes[i]); std::set<MVertex *>::iterator itN = myNodes.find(bcNodes[i]);
if (itN != myNodes.end()){ if (itN != myNodes.end()){
std::map<MVertex*, int>::iterator itm = _mapV.find(bcNodes[i]); std::map<MVertex*, int>::iterator itm = _mapV.find(bcNodes[i]);
...@@ -954,7 +954,7 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper, ...@@ -954,7 +954,7 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
F.scale(0.0); F.scale(0.0);
//UNIT CIRCLE //UNIT CIRCLE
for (int i=0; i < bcNodes.size(); i++){ for (unsigned int i = 0; i < bcNodes.size(); i++){
std::set<MVertex *>::iterator itN = myNodes.find(bcNodes[i]); std::set<MVertex *>::iterator itN = myNodes.find(bcNodes[i]);
if (itN != myNodes.end()){ if (itN != myNodes.end()){
std::map<MVertex*, int>::iterator itm = _mapV.find(bcNodes[i]); std::map<MVertex*, int>::iterator itm = _mapV.find(bcNodes[i]);
...@@ -1112,159 +1112,3 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval, ...@@ -1112,159 +1112,3 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
lastDXDV = dXdv; lastDXDV = dXdv;
return true; return true;
} }
// bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
// double &XX, double &YY, double &ZZ,
// SVector3 &dXdu, SVector3& dXdv, int num_neighbours){
// //Finds the U,V,S (in the 0-level set) that are the 'num_neighbours' closest to u_eval and v_eval.
// //Thus in total, we're working with '3*num_neighbours' nodes
// //Say that the vector 'index' gives us the indices of the closest points
// bool converged = true;
// num_neighbours = 30;
// #if defined (HAVE_ANN)
// double uvw[3] = { u_eval, v_eval, 0.0 };
// ANNidxArray index = new ANNidx[num_neighbours];
// ANNdistArray dist = new ANNdist[num_neighbours];
// UVkdtree->annkSearch(uvw, num_neighbours, index, dist);
// int N_nbr = 5;//num_neighbours;
// fullMatrix<double> ux(N_nbr,3), uy(N_nbr,3), uz(N_nbr,3), u_temp(N_nbr,3);
// fullMatrix<double> nodes_eval(N_nbr,3),temp_nodes_eval(1,3),temp_u_temp(1,3);
// fullMatrix<double> u_vec(3*num_neighbours,3);
// fullMatrix<double> xyz_local(3*num_neighbours,3);
// // u_vec = [u v s] : These are the u,v,s that are on the surface *and* outside of it also!
// for (int i = 0; i < num_neighbours; i++){
// u_vec(i,0) = UV(index[i],0);
// u_vec(i,1) = UV(index[i],1);
// u_vec(i,2) = 0.0;
// u_vec(i+num_neighbours,0) = UV(index[i]+nbNodes,0);
// u_vec(i+num_neighbours,1) = UV(index[i]+nbNodes,1);
// u_vec(i+num_neighbours,2) = surfInterp(index[i]+nbNodes,0)*deltaUV;
// u_vec(i+2*num_neighbours,0) = UV(index[i]+2*nbNodes,0);
// u_vec(i+2*num_neighbours,1) = UV(index[i]+2*nbNodes,1);
// u_vec(i+2*num_neighbours,2) = surfInterp(index[i]+2*nbNodes,0)*deltaUV;
// xyz_local(i,0) = extendedX(index[i],0);
// xyz_local(i,1) = extendedX(index[i],1);
// xyz_local(i,2) = extendedX(index[i],2);
// xyz_local(i+num_neighbours,0) = extendedX(index[i]+nbNodes,0);
// xyz_local(i+num_neighbours,1) = extendedX(index[i]+nbNodes,1);
// xyz_local(i+num_neighbours,2) = extendedX(index[i]+nbNodes,2);
// xyz_local(i+2*num_neighbours,0) = extendedX(index[i]+2*nbNodes,0);
// xyz_local(i+2*num_neighbours,1) = extendedX(index[i]+2*nbNodes,1);
// xyz_local(i+2*num_neighbours,2) = extendedX(index[i]+2*nbNodes,2);
// }
// // u_vec_eval = [u_eval v_eval s_eval]
// fullMatrix<double> u_vec_eval(1, 3),nodes_vec_eval(1,3),u_test_vec_eval(1,3);
// u_vec_eval(0,0) = u_eval;
// u_vec_eval(0,1) = v_eval;
// u_vec_eval(0,2) = 0.0;
// //we will use a local interpolation to find the corresponding XYZ point to (u_eval,v_eval).
// _inUV = 1;
// evalRbfDer(0, u_vec, u_vec_eval,xyz_local, temp_nodes_eval);
// nodes_eval(0,0) = temp_nodes_eval(0,0);
// nodes_eval(0,1) = temp_nodes_eval(0,1);
// nodes_eval(0,2) = temp_nodes_eval(0,2);
// _inUV= 0;
// evalRbfDer(0,xyz_local,nodes_eval,u_vec,temp_u_temp);
// u_temp(0,0) = temp_u_temp(0,0);
// u_temp(0,1) = temp_u_temp(0,1);
// u_temp(0,2) = temp_u_temp(0,2);
// //we use the closest point
// for (int i_nbr = 1; i_nbr < N_nbr; i_nbr++){
// nodes_eval(i_nbr,0) = extendedX(index[i_nbr-1],0);
// nodes_eval(i_nbr,1) = extendedX(index[i_nbr-1],1);
// nodes_eval(i_nbr,2) = extendedX(index[i_nbr-1],2);
// u_temp(i_nbr,0) = UV(index[i_nbr-1],0);
// u_temp(i_nbr,1) = UV(index[i_nbr-1],1);
// u_temp(i_nbr,2) = 0.0;
// }
// delete [] index;
// delete [] dist;
// int incr = 0;
// double norm = 0.0;
// double norm_temp = 0.0;
// double dist_test = 0.0;
// fullMatrix<double> Jac(3,3),best_Jac(3,3);
// fullMatrix<double> normDist(N_nbr,1);
// fullMatrix<double> nodes_eval_temp(N_nbr,3);
// double norm_treshhold = 8.0;
// std::vector<fullMatrix<double> >JacV;
// while(norm < norm_treshhold && incr < 10){
// // Find the entries of the m Jacobians
// evalRbfDer(1,xyz_local,nodes_eval,u_vec,ux);
// evalRbfDer(2,xyz_local,nodes_eval,u_vec,uy);
// evalRbfDer(3,xyz_local,nodes_eval,u_vec,uz);
// for (int i_nbr = 0; i_nbr < N_nbr; i_nbr++){
// // Jacobian of the UVS coordinate system w.r.t. XYZ
// for (int k = 0; k < 3;k++){
// Jac(k,0) = ux(i_nbr,k);
// Jac(k,1) = uy(i_nbr,k);
// Jac(k,2) = uz(i_nbr,k);
// }
// Jac.invertInPlace();
// JacV.push_back(Jac);
// for (int j = 0; j< 3;j++)
// nodes_eval(i_nbr,j) = nodes_eval(i_nbr,j)
// + Jac(j,0)*( u_vec_eval(0,0) - u_temp(i_nbr,0))
// + Jac(j,1)*( u_vec_eval(0,1) - u_temp(i_nbr,1))
// + Jac(j,2)*( 0.0 - u_temp(i_nbr,2));
// }
// // Find the corresponding u, v, s
// evalRbfDer(0,xyz_local,nodes_eval,u_vec,u_temp);
// norm = 0;
// for (int i_nbr = 0; i_nbr < N_nbr; i_nbr++){
// normDist(i_nbr,0) = sqrt((u_temp(i_nbr,0)-u_vec_eval(0,0))*(u_temp(i_nbr,0)-u_vec_eval(0,0))+(u_temp(i_nbr,1)-u_vec_eval(0,1))*(u_temp(i_nbr,1)-u_vec_eval(0,1))+(u_temp(i_nbr,2)*u_temp(i_nbr,2)));
// norm_temp = -(log10(normDist(i_nbr,0)));
// if (norm_temp>norm_treshhold){
// norm = norm_temp;
// nodes_vec_eval(0,0)=nodes_eval(i_nbr,0);
// nodes_vec_eval(0,1)=nodes_eval(i_nbr,1);
// nodes_vec_eval(0,2)=nodes_eval(i_nbr,2);
// u_test_vec_eval(0,0)=u_temp(i_nbr,0);
// u_test_vec_eval(0,1)=u_temp(i_nbr,1);
// u_test_vec_eval(0,2)=u_temp(i_nbr,2);
// best_Jac = JacV[i_nbr];
// dist_test = normDist(i_nbr,0);
// i_nbr=N_nbr; //To get out of the for loof and thus also to get out of the while loop
// }
// }
// incr++;
// }
// if (norm < norm_treshhold ){
// printf("Newton not converged (norm=%g, dist=%g) for uv=(%g,%g) UVS=(%g,%g,%g) NB=%d \n", norm, dist_test,u_eval, v_eval, u_test_vec_eval(0,0), u_test_vec_eval(0,1), u_test_vec_eval(0,2), num_neighbours);
// if (num_neighbours+5 < nbNodes)
// return UVStoXYZ(u_eval, v_eval, XX, YY,ZZ, dXdu, dXdv, num_neighbours+5);
// else converged = false;
// }
// XX = nodes_vec_eval(0,0)*sBox;
// YY = nodes_vec_eval(0,1)*sBox;
// ZZ = nodes_vec_eval(0,2)*sBox;
// dXdu = SVector3(best_Jac(0,0)*sBox, best_Jac(1,0)*sBox, best_Jac(2,0)*sBox);
// dXdv = SVector3(best_Jac(0,1)*sBox, best_Jac(1,1)*sBox, best_Jac(2,1)*sBox);
// return converged;
// #endif
// }
Geo/gmshLevelset.cpp
+ 258
− 153
View file @ e02840df
...@@ -31,6 +31,7 @@ void insertActiveCells(double x, double y, double z, double rmax, ...@@ -31,6 +31,7 @@ void insertActiveCells(double x, double y, double z, double rmax,
for (int k = k1; k <= k2; k++) for (int k = k1; k <= k2; k++)
box.insertActiveCell(box.getCellIndex(i, j, k)); box.insertActiveCell(box.getCellIndex(i, j, k));
} }
void fillPointCloud(GEdge *ge, double sampling, std::vector<SPoint3> &points) void fillPointCloud(GEdge *ge, double sampling, std::vector<SPoint3> &points)
{ {
Range<double> t_bounds = ge->parBounds(0); Range<double> t_bounds = ge->parBounds(0);
...@@ -108,8 +109,6 @@ void removeParentCellsWithChildren(cartesianBox<double> *box) ...@@ -108,8 +109,6 @@ void removeParentCellsWithChildren(cartesianBox<double> *box)
void computeLevelset(GModel *gm, cartesianBox<double> &box) void computeLevelset(GModel *gm, cartesianBox<double> &box)
{ {
// tolerance for desambiguation
const double tol = box.getLC() * 1.e-12;
std::vector<SPoint3> nodes; std::vector<SPoint3> nodes;
std::vector<int> indices; std::vector<int> indices;
for (cartesianBox<double>::valIter it = box.nodalValuesBegin(); for (cartesianBox<double>::valIter it = box.nodalValuesBegin();
...@@ -160,21 +159,24 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box) ...@@ -160,21 +159,24 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
if(box.getChildBox()) computeLevelset(gm, *box.getChildBox()); if(box.getChildBox()) computeLevelset(gm, *box.getChildBox());
} }
inline double det3(double d11, double d12, double d13,
double d21, double d22, double d23,
double d31, double d32, double d33)
inline double det3(double d11, double d12, double d13, double d21, double d22, double d23, double d31, double d32, double d33) { {
return d11 * (d22 * d33 - d23 * d32) - d21 * (d12 * d33 - d13 * d32) + d31 * (d12 * d23 - d13 * d22); return d11 * (d22 * d33 - d23 * d32) - d21 * (d12 * d33 - d13 * d32) +
d31 * (d12 * d23 - d13 * d22);
} }
inline void norm(const double *vec, double *norm) { inline void norm(const double *vec, double *norm)
{
double n = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); double n = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
norm[0] = vec[0] / n; norm[0] = vec[0] / n;
norm[1] = vec[1] / n; norm[1] = vec[1] / n;
norm[2] = vec[2] / n; norm[2] = vec[2] / n;
} }
inline void cross(const double *pt0, const double *pt1, const double *pt2, double *cross) { inline void cross(const double *pt0, const double *pt1, const double *pt2, double *cross)
{
double v1[3] = {pt1[0] - pt0[0], pt1[1] - pt0[1], pt1[2] - pt0[2]}; double v1[3] = {pt1[0] - pt0[0], pt1[1] - pt0[1], pt1[2] - pt0[2]};
double v2[3] = {pt2[0] - pt0[0], pt2[1] - pt0[1], pt2[2] - pt0[2]}; double v2[3] = {pt2[0] - pt0[0], pt2[1] - pt0[1], pt2[2] - pt0[2]};
cross[0] = v1[1] * v2[2] - v2[1] * v1[2]; cross[0] = v1[1] * v2[2] - v2[1] * v1[2];
...@@ -182,7 +184,9 @@ inline void cross(const double *pt0, const double *pt1, const double *pt2, doubl ...@@ -182,7 +184,9 @@ inline void cross(const double *pt0, const double *pt1, const double *pt2, doubl
cross[2] = v1[0] * v2[1] - v2[0] * v1[1]; cross[2] = v1[0] * v2[1] - v2[0] * v1[1];
} }
inline bool isPlanar(const double *pt1, const double *pt2, const double *pt3, const double *pt4) { inline bool isPlanar(const double *pt1, const double *pt2, const double *pt3,
const double *pt4)
{
double c1[3]; cross(pt1, pt2, pt3, c1); double c1[3]; cross(pt1, pt2, pt3, c1);
double c2[3]; cross(pt1, pt2, pt4, c2); double c2[3]; cross(pt1, pt2, pt4, c2);
double n1[3]; norm(c1, n1); double n1[3]; norm(c1, n1);
...@@ -190,8 +194,9 @@ inline bool isPlanar(const double *pt1, const double *pt2, const double *pt3, co ...@@ -190,8 +194,9 @@ inline bool isPlanar(const double *pt1, const double *pt2, const double *pt3, co
return (n1[0] == n2[0] && n1[1] == n2[1] && n1[2] == n2[2]); return (n1[0] == n2[0] && n1[1] == n2[1] && n1[2] == n2[2]);
} }
inline double evalRadialFnDer (int p, int index, double dx, double dy, double dz, double ep){ inline double evalRadialFnDer (int p, int index, double dx, double dy, double dz,
double ep)
{
double r2 = dx*dx+dy*dy+dz*dz; //r^2 double r2 = dx*dx+dy*dy+dz*dz; //r^2
switch (index) { switch (index) {
case 0 : // GA case 0 : // GA
...@@ -209,20 +214,24 @@ inline double evalRadialFnDer (int p, int index, double dx, double dy, double dz ...@@ -209,20 +214,24 @@ inline double evalRadialFnDer (int p, int index, double dx, double dy, double dz
case 3: return ep*ep*dz/sqrt(1+ep*ep*r2); case 3: return ep*ep*dz/sqrt(1+ep*ep*r2);
} }
} }
return 0.;
} }
inline void printNodes(fullMatrix<double> &myNodes, fullMatrix<double> &surf){ inline void printNodes(fullMatrix<double> &myNodes, fullMatrix<double> &surf)
{
FILE * xyz = fopen("myNodes.pos","w"); FILE * xyz = fopen("myNodes.pos","w");
fprintf(xyz,"View \"\"{\n"); fprintf(xyz,"View \"\"{\n");
for(int itv = 1; itv != myNodes.size1(); ++itv){ for(int itv = 1; itv != myNodes.size1(); ++itv){
fprintf(xyz,"SP(%g,%g,%g){%g};\n", myNodes(itv,0),myNodes(itv,1),myNodes(itv,2),surf(itv,0)); fprintf(xyz,"SP(%g,%g,%g){%g};\n", myNodes(itv,0), myNodes(itv,1),
myNodes(itv,2), surf(itv,0));
} }
fprintf(xyz,"};\n"); fprintf(xyz,"};\n");
fclose(xyz); fclose(xyz);
} }
// extrude a list of the primitive levelsets with a "Level-order traversal sequence" // extrude a list of the primitive levelsets with a "Level-order traversal sequence"
void gLevelset::getPrimitives(std::vector<gLevelset *> &gLsPrimitives) { void gLevelset::getPrimitives(std::vector<gLevelset *> &gLsPrimitives)
{
std::queue<gLevelset *> Q; std::queue<gLevelset *> Q;
Q.push(this); Q.push(this);
while(!Q.empty()){ while(!Q.empty()){
...@@ -238,8 +247,10 @@ void gLevelset::getPrimitives(std::vector<gLevelset *> &gLsPrimitives) { ...@@ -238,8 +247,10 @@ void gLevelset::getPrimitives(std::vector<gLevelset *> &gLsPrimitives) {
} }
} }
} }
// extrude a list of the primitive levelsets with a "post-order traversal sequence" // extrude a list of the primitive levelsets with a "post-order traversal sequence"
void gLevelset::getPrimitivesPO(std::vector<gLevelset *> &gLsPrimitives) { void gLevelset::getPrimitivesPO(std::vector<gLevelset *> &gLsPrimitives)
{
std::stack<gLevelset *> S; std::stack<gLevelset *> S;
std::stack<gLevelset *> Sc; // levelset checked std::stack<gLevelset *> Sc; // levelset checked
S.push(this); S.push(this);
...@@ -268,7 +279,8 @@ void gLevelset::getPrimitivesPO(std::vector<gLevelset *> &gLsPrimitives) { ...@@ -268,7 +279,8 @@ void gLevelset::getPrimitivesPO(std::vector<gLevelset *> &gLsPrimitives) {
} }
// return a list with the levelsets in a "Reverse Polish Notation" // return a list with the levelsets in a "Reverse Polish Notation"
void gLevelset::getRPN(std::vector<gLevelset *> &gLsRPN) { void gLevelset::getRPN(std::vector<gLevelset *> &gLsRPN)
{
std::stack<gLevelset *> S; std::stack<gLevelset *> S;
std::stack<gLevelset *> Sc; // levelset checked std::stack<gLevelset *> Sc; // levelset checked
S.push(this); S.push(this);
...@@ -302,25 +314,37 @@ gLevelset::gLevelset(const gLevelset &lv) ...@@ -302,25 +314,37 @@ gLevelset::gLevelset(const gLevelset &lv)
tag_ = lv.tag_; tag_ = lv.tag_;
} }
gLevelsetPlane::gLevelsetPlane(const std::vector<double> &pt, const std::vector<double> &norm, int tag) : gLevelsetPrimitive(tag) { gLevelsetPlane::gLevelsetPlane(const std::vector<double> &pt,
const std::vector<double> &norm, int tag)
: gLevelsetPrimitive(tag)
{
a = norm[0]; a = norm[0];
b = norm[1]; b = norm[1];
c = norm[2]; c = norm[2];
d = -a * pt[0] - b * pt[1] - c * pt[2]; d = -a * pt[0] - b * pt[1] - c * pt[2];
} }
gLevelsetPlane::gLevelsetPlane(const double * pt, const double *norm, int tag) : gLevelsetPrimitive(tag) {
gLevelsetPlane::gLevelsetPlane(const double * pt, const double *norm, int tag)
: gLevelsetPrimitive(tag)
{
a = norm[0]; a = norm[0];
b = norm[1]; b = norm[1];
c = norm[2]; c = norm[2];
d = -a * pt[0] - b * pt[1] - c * pt[2]; d = -a * pt[0] - b * pt[1] - c * pt[2];
} }
gLevelsetPlane::gLevelsetPlane(const double * pt1, const double *pt2, const double *pt3, int tag) : gLevelsetPrimitive(tag) { gLevelsetPlane::gLevelsetPlane(const double * pt1, const double *pt2,
const double *pt3, int tag)
: gLevelsetPrimitive(tag)
{
a = det3(1., pt1[1], pt1[2], 1., pt2[1], pt2[2], 1., pt3[1], pt3[2]); a = det3(1., pt1[1], pt1[2], 1., pt2[1], pt2[2], 1., pt3[1], pt3[2]);
b = det3(pt1[0], 1., pt1[2], pt2[0], 1., pt2[2], pt3[0], 1., pt3[2]); b = det3(pt1[0], 1., pt1[2], pt2[0], 1., pt2[2], pt3[0], 1., pt3[2]);
c = det3(pt1[0], pt1[1], 1., pt2[0], pt2[1], 1., pt3[0], pt3[1], 1.); c = det3(pt1[0], pt1[1], 1., pt2[0], pt2[1], 1., pt3[0], pt3[1], 1.);
d = -det3(pt1[0], pt1[1], pt1[2], pt2[0], pt2[1], pt2[2], pt3[0], pt3[1], pt3[2]); d = -det3(pt1[0], pt1[1], pt1[2], pt2[0], pt2[1], pt2[2], pt3[0], pt3[1], pt3[2]);
} }
gLevelsetPlane::gLevelsetPlane(const gLevelsetPlane &lv) : gLevelsetPrimitive(lv) {
gLevelsetPlane::gLevelsetPlane(const gLevelsetPlane &lv)
: gLevelsetPrimitive(lv)
{
a = lv.a; a = lv.a;
b = lv.b; b = lv.b;
c = lv.c; c = lv.c;
...@@ -330,7 +354,8 @@ gLevelsetPlane::gLevelsetPlane(const gLevelsetPlane &lv) : gLevelsetPrimitive(lv ...@@ -330,7 +354,8 @@ gLevelsetPlane::gLevelsetPlane(const gLevelsetPlane &lv) : gLevelsetPrimitive(lv
//level set defined by points (RBF interpolation) //level set defined by points (RBF interpolation)
fullMatrix<double> gLevelsetPoints::generateRbfMat(int p, int index, fullMatrix<double> gLevelsetPoints::generateRbfMat(int p, int index,
const fullMatrix<double> &nodes1, const fullMatrix<double> &nodes1,
const fullMatrix<double> &nodes2) const { const fullMatrix<double> &nodes2) const
{
int m = nodes2.size1(); int m = nodes2.size1();
int n = nodes1.size1(); int n = nodes1.size1();
fullMatrix<double> rbfMat(m,n); fullMatrix<double> rbfMat(m,n);
...@@ -347,12 +372,13 @@ fullMatrix<double> gLevelsetPoints::generateRbfMat(int p, int index, ...@@ -347,12 +372,13 @@ fullMatrix<double> gLevelsetPoints::generateRbfMat(int p, int index,
return rbfMat; return rbfMat;
} }
void gLevelsetPoints::RbfOp(int p, int index, void gLevelsetPoints::RbfOp(int p, int index,
const fullMatrix<double> &cntrs, const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes, const fullMatrix<double> &nodes,
fullMatrix<double> &D, fullMatrix<double> &D,
bool isLocal) const { bool isLocal) const
{
fullMatrix<double> rbfMatB = generateRbfMat(p,index, cntrs,nodes); fullMatrix<double> rbfMatB = generateRbfMat(p,index, cntrs,nodes);
//printf("size=%d %d \n", rbfMatB.size1(), rbfMatB.size2()); //printf("size=%d %d \n", rbfMatB.size1(), rbfMatB.size2());
//rbfMatB.print("MatB"); //rbfMatB.print("MatB");
...@@ -373,32 +399,30 @@ void gLevelsetPoints::RbfOp(int p, int index, ...@@ -373,32 +399,30 @@ void gLevelsetPoints::RbfOp(int p, int index,
} }
void gLevelsetPoints::evalRbfDer(int p, int index, void gLevelsetPoints::evalRbfDer(int p, int index,
const fullMatrix<double> &cntrs, const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes, const fullMatrix<double> &nodes,
const fullMatrix<double> &fValues, const fullMatrix<double> &fValues,
fullMatrix<double> &fApprox, bool isLocal) const { fullMatrix<double> &fApprox, bool isLocal) const
{
fApprox.resize(nodes.size1(),fValues.size2()); fApprox.resize(nodes.size1(),fValues.size2());
fullMatrix<double> D; fullMatrix<double> D;
RbfOp(p,index, cntrs,nodes,D,isLocal); RbfOp(p,index, cntrs,nodes,D,isLocal);
//D.print("D"); //D.print("D");
fApprox.gemm(D,fValues, 1.0, 0.0); fApprox.gemm(D,fValues, 1.0, 0.0);
} }
void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs, void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs,
fullMatrix<double> &level_set_nodes, fullMatrix<double> &level_set_nodes,
fullMatrix<double> &level_set_funvals){ fullMatrix<double> &level_set_funvals)
{
int numNodes = cntrs.size1(); int numNodes = cntrs.size1();
int nTot = 3*numNodes; int nTot = 3*numNodes;
double normFactor; double normFactor;
level_set_nodes.resize(nTot,3); level_set_nodes.resize(nTot,3);
level_set_funvals.resize(nTot,1); level_set_funvals.resize(nTot,1);
fullMatrix<double> ONES(numNodes+1,1),sx(numNodes,1),sy(numNodes,1),sz(numNodes,1),norms(numNodes,3), cntrsPlus(numNodes+1,3); fullMatrix<double> ONES(numNodes+1,1),sx(numNodes,1),sy(numNodes,1);
fullMatrix<double> sz(numNodes,1),norms(numNodes,3), cntrsPlus(numNodes+1,3);
//Computes the normal vectors to the surface at each node //Computes the normal vectors to the surface at each node
double dist_min = 1.e6; double dist_min = 1.e6;
...@@ -445,8 +469,9 @@ void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs, ...@@ -445,8 +469,9 @@ void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs,
} }
} }
gLevelsetPoints::gLevelsetPoints(fullMatrix<double> &centers, int tag)
gLevelsetPoints::gLevelsetPoints(fullMatrix<double> &centers, int tag) : gLevelsetPrimitive(tag) { : gLevelsetPrimitive(tag)
{
int nbNodes = 3*centers.size1(); int nbNodes = 3*centers.size1();
setup_level_set(centers, points, surf); setup_level_set(centers, points, surf);
...@@ -459,15 +484,16 @@ gLevelsetPoints::gLevelsetPoints(fullMatrix<double> &centers, int tag) : gLevels ...@@ -459,15 +484,16 @@ gLevelsetPoints::gLevelsetPoints(fullMatrix<double> &centers, int tag) : gLevels
matAInv.invertInPlace(); matAInv.invertInPlace();
//printf("End init levelset points %d \n", points.size1()); //printf("End init levelset points %d \n", points.size1());
} }
gLevelsetPoints::gLevelsetPoints(const gLevelsetPoints &lv) : gLevelsetPrimitive(lv) { gLevelsetPoints::gLevelsetPoints(const gLevelsetPoints &lv)
: gLevelsetPrimitive(lv)
{
points = lv.points; points = lv.points;
} }
double gLevelsetPoints::operator()(const double x, const double y, const double z) const{ double gLevelsetPoints::operator()(const double x, const double y, const double z) const
{
if(mapP.empty()) printf("Levelset Points : call computeLS() before calling operator()\n"); if(mapP.empty()) printf("Levelset Points : call computeLS() before calling operator()\n");
SPoint3 sp(x,y,z); SPoint3 sp(x,y,z);
...@@ -483,19 +509,18 @@ double gLevelsetPoints::operator()(const double x, const double y, const double ...@@ -483,19 +509,18 @@ double gLevelsetPoints::operator()(const double x, const double y, const double
// xyz_eval(0,2) = z; // xyz_eval(0,2) = z;
// evalRbfDer(0, 1, points, xyz_eval, surf, surf_eval); // evalRbfDer(0, 1, points, xyz_eval, surf, surf_eval);
// return surf_eval(0,0); // return surf_eval(0,0);
} }
void gLevelsetPoints::computeLS(std::vector<MVertex*> &vert){ void gLevelsetPoints::computeLS(std::vector<MVertex*> &vert)
{
fullMatrix<double> xyz_eval(vert.size(), 3), surf_eval(vert.size(), 1); fullMatrix<double> xyz_eval(vert.size(), 3), surf_eval(vert.size(), 1);
for (int i = 0; i< vert.size(); i++){ for (unsigned int i = 0; i < vert.size(); i++){
xyz_eval(i,0) = vert[i]->x(); xyz_eval(i,0) = vert[i]->x();
xyz_eval(i,1) = vert[i]->y(); xyz_eval(i,1) = vert[i]->y();
xyz_eval(i,2) = vert[i]->z(); xyz_eval(i,2) = vert[i]->z();
} }
evalRbfDer(0, 1, points, xyz_eval, surf, surf_eval); evalRbfDer(0, 1, points, xyz_eval, surf, surf_eval);
for (int i = 0; i< vert.size(); i++){ for (unsigned int i = 0; i < vert.size(); i++){
mapP[SPoint3(vert[i]->x(), vert[i]->y(),vert[i]->z())] = surf_eval(i,0); mapP[SPoint3(vert[i]->x(), vert[i]->y(),vert[i]->z())] = surf_eval(i,0);
} }
} }
...@@ -511,7 +536,9 @@ void gLevelsetPoints::computeLS(std::vector<MVertex*> &vert){ ...@@ -511,7 +536,9 @@ void gLevelsetPoints::computeLS(std::vector<MVertex*> &vert){
x^T A x + [b^T - 2 A t] x + [c - b^T t + t^T A t ] = 0 x^T A x + [b^T - 2 A t] x + [c - b^T t + t^T A t ] = 0
*/ */
gLevelsetQuadric::gLevelsetQuadric(const gLevelsetQuadric &lv) : gLevelsetPrimitive(lv){ gLevelsetQuadric::gLevelsetQuadric(const gLevelsetQuadric &lv)
: gLevelsetPrimitive(lv)
{
for(int i = 0; i < 3; i++){ for(int i = 0; i < 3; i++){
B[i] = lv.B[i]; B[i] = lv.B[i];
for(int j = 0; j < 3; j++) for(int j = 0; j < 3; j++)
...@@ -519,7 +546,9 @@ gLevelsetQuadric::gLevelsetQuadric(const gLevelsetQuadric &lv) : gLevelsetPrimit ...@@ -519,7 +546,9 @@ gLevelsetQuadric::gLevelsetQuadric(const gLevelsetQuadric &lv) : gLevelsetPrimit
} }
C = lv.C; C = lv.C;
} }
void gLevelsetQuadric::Ax(const double x[3], double res[3], double fact){
void gLevelsetQuadric::Ax(const double x[3], double res[3], double fact)
{
for(int i = 0; i < 3; i++){ for(int i = 0; i < 3; i++){
res[i] = 0.; res[i] = 0.;
for(int j = 0; j < 3; j++){ for(int j = 0; j < 3; j++){
...@@ -528,12 +557,15 @@ void gLevelsetQuadric::Ax(const double x[3], double res[3], double fact){ ...@@ -528,12 +557,15 @@ void gLevelsetQuadric::Ax(const double x[3], double res[3], double fact){
} }
} }
void gLevelsetQuadric::xAx(const double x[3], double &res, double fact){ void gLevelsetQuadric::xAx(const double x[3], double &res, double fact)
res= fact * (A[0][0] * x[0] * x[0] + A[1][1] * x[1] * x[1] + A[2][2] * x[2] * x[2] {
+ A[1][0] * x[1] * x[0] * 2. + A[2][0] * x[2] * x[0] * 2. + A[1][2] * x[1] * x[2] * 2.); res= fact * (A[0][0] * x[0] * x[0] + A[1][1] * x[1] * x[1] +
A[2][2] * x[2] * x[2] + A[1][0] * x[1] * x[0] * 2. +
A[2][0] * x[2] * x[0] * 2. + A[1][2] * x[1] * x[2] * 2.);
} }
void gLevelsetQuadric::translate(const double transl[3]){ void gLevelsetQuadric::translate(const double transl[3])
{
double b; double b;
xAx(transl, b, 1.0); xAx(transl, b, 1.0);
C += (-B[0] * transl[0] - B[1] * transl[1] - B[2] * transl[2] + b); C += (-B[0] * transl[0] - B[1] * transl[1] - B[2] * transl[2] + b);
...@@ -545,7 +577,8 @@ void gLevelsetQuadric::translate(const double transl[3]){ ...@@ -545,7 +577,8 @@ void gLevelsetQuadric::translate(const double transl[3]){
B[2] += -x[2]; B[2] += -x[2];
} }
void gLevelsetQuadric::rotate(const double rotate[3][3]){ void gLevelsetQuadric::rotate(const double rotate[3][3])
{
double a11 = 0., a12 = 0., a13 = 0., a22 = 0., a23 = 0., a33 = 0., b1 = 0., b2 = 0., b3 = 0.; double a11 = 0., a12 = 0., a13 = 0., a22 = 0., a23 = 0., a33 = 0., b1 = 0., b2 = 0., b3 = 0.;
for(int i = 0; i < 3; i++){ for(int i = 0; i < 3; i++){
b1 += B[i] * rotate[i][0]; b1 += B[i] * rotate[i][0];
...@@ -572,7 +605,8 @@ void gLevelsetQuadric::rotate(const double rotate[3][3]){ ...@@ -572,7 +605,8 @@ void gLevelsetQuadric::rotate(const double rotate[3][3]){
} }
// computes the rotation matrix of the rotation of the vector (0,0,1) to dir // computes the rotation matrix of the rotation of the vector (0,0,1) to dir
void gLevelsetQuadric::computeRotationMatrix(const double dir[3], double t[3][3]){ void gLevelsetQuadric::computeRotationMatrix(const double dir[3], double t[3][3])
{
double norm = sqrt(dir[0] * dir[0] + dir[1] * dir[1]); double norm = sqrt(dir[0] * dir[0] + dir[1] * dir[1]);
double n[3] = {1., 0., 0.}; double n[3] = {1., 0., 0.};
double ct = 1., st = 0.; double ct = 1., st = 0.;
...@@ -593,10 +627,16 @@ void gLevelsetQuadric::computeRotationMatrix(const double dir[3], double t[3][3] ...@@ -593,10 +627,16 @@ void gLevelsetQuadric::computeRotationMatrix(const double dir[3], double t[3][3]
t[2][2] = n[2] * n[2] * (1. - ct) + ct; t[2][2] = n[2] * n[2] * (1. - ct) + ct;
} }
void gLevelsetQuadric::computeRotationMatrix(const double dir1[3], const double dir2[3], double t[3][3]){ void gLevelsetQuadric::computeRotationMatrix(const double dir1[3],
double norm = sqrt((dir1[1] * dir2[2] - dir1[2] * dir2[1]) * (dir1[1] * dir2[2] - dir1[2] * dir2[1]) const double dir2[3],
+ (dir1[2] * dir2[0] - dir1[0] * dir2[2]) * (dir1[2] * dir2[0] - dir1[0] * dir2[2]) double t[3][3])
+ (dir1[0] * dir2[1] - dir1[1] * dir2[0]) * (dir1[0] * dir2[1] - dir1[1] * dir2[0])); {
double norm = sqrt((dir1[1] * dir2[2] - dir1[2] * dir2[1]) *
(dir1[1] * dir2[2] - dir1[2] * dir2[1])
+ (dir1[2] * dir2[0] - dir1[0] * dir2[2]) *
(dir1[2] * dir2[0] - dir1[0] * dir2[2])
+ (dir1[0] * dir2[1] - dir1[1] * dir2[0]) *
(dir1[0] * dir2[1] - dir1[1] * dir2[0]));
double n[3] = {1., 0., 0.}; double n[3] = {1., 0., 0.};
double ct = 1., st = 0.; double ct = 1., st = 0.;
if(norm != 0.) { if(norm != 0.) {
...@@ -618,7 +658,8 @@ void gLevelsetQuadric::computeRotationMatrix(const double dir1[3], const double ...@@ -618,7 +658,8 @@ void gLevelsetQuadric::computeRotationMatrix(const double dir1[3], const double
t[2][2] = n[2] * n[2] * (1. - ct) + ct; t[2][2] = n[2] * n[2] * (1. - ct) + ct;
} }
void gLevelsetQuadric::init(){ void gLevelsetQuadric::init()
{
for(int i = 0; i < 3; i++){ for(int i = 0; i < 3; i++){
B[i] = 0.; B[i] = 0.;
for(int j = 0; j < 3; j++) A[i][j] = 0.; for(int j = 0; j < 3; j++) A[i][j] = 0.;
...@@ -626,12 +667,16 @@ void gLevelsetQuadric::init(){ ...@@ -626,12 +667,16 @@ void gLevelsetQuadric::init(){
C = 0.; C = 0.;
} }
double gLevelsetQuadric::operator()(const double x, const double y, const double z) const{ double gLevelsetQuadric::operator()(const double x, const double y, const double z) const
{
return(A[0][0] * x * x + 2. * A[0][1] * x * y + 2. * A[0][2] * x * z + A[1][1] * y * y return(A[0][0] * x * x + 2. * A[0][1] * x * y + 2. * A[0][2] * x * z + A[1][1] * y * y
+ 2. * A[1][2] * y * z + A[2][2] * z * z + B[0] * x + B[1] * y + B[2] * z + C); + 2. * A[1][2] * y * z + A[2][2] * z * z + B[0] * x + B[1] * y + B[2] * z + C);
} }
gLevelsetShamrock::gLevelsetShamrock(double _xmid, double _ymid, double _zmid, double _a, double _b, int _c, int tag) : gLevelsetPrimitive(tag), xmid(_xmid), ymid(_ymid), zmid(_zmid), a(_a), b(_b), c(_c) { gLevelsetShamrock::gLevelsetShamrock(double _xmid, double _ymid, double _zmid,
double _a, double _b, int _c, int tag)
: gLevelsetPrimitive(tag), xmid(_xmid), ymid(_ymid), zmid(_zmid), a(_a), b(_b), c(_c)
{
// creating the iso-zero // creating the iso-zero
double angle = 0.; double angle = 0.;
double r; double r;
...@@ -644,7 +689,8 @@ gLevelsetShamrock::gLevelsetShamrock(double _xmid, double _ymid, double _zmid, d ...@@ -644,7 +689,8 @@ gLevelsetShamrock::gLevelsetShamrock(double _xmid, double _ymid, double _zmid, d
} }
} }
double gLevelsetShamrock::operator() (const double x, const double y, const double z) const { double gLevelsetShamrock::operator() (const double x, const double y, const double z) const
{
// computing distance to pre-defined (sampled) iso-zero // computing distance to pre-defined (sampled) iso-zero
double dx,dy,xi,yi,d; double dx,dy,xi,yi,d;
dx = x-iso_x[0]; dx = x-iso_x[0];
...@@ -687,7 +733,10 @@ double gLevelsetShamrock::operator() (const double x, const double y, const doub ...@@ -687,7 +733,10 @@ double gLevelsetShamrock::operator() (const double x, const double y, const doub
return (sign*distance); return (sign*distance);
} }
gLevelsetPopcorn::gLevelsetPopcorn(double _xc, double _yc, double _zc, double _r0, double _A, double _sigma, int tag) : gLevelsetPrimitive(tag) { gLevelsetPopcorn::gLevelsetPopcorn(double _xc, double _yc, double _zc, double _r0,
double _A, double _sigma, int tag)
: gLevelsetPrimitive(tag)
{
A = _A; A = _A;
sigma = _sigma; sigma = _sigma;
r0 = _r0; r0 = _r0;
...@@ -696,7 +745,8 @@ gLevelsetPopcorn::gLevelsetPopcorn(double _xc, double _yc, double _zc, double _r ...@@ -696,7 +745,8 @@ gLevelsetPopcorn::gLevelsetPopcorn(double _xc, double _yc, double _zc, double _r
zc = _zc; zc = _zc;
} }
double gLevelsetPopcorn::operator() (const double x, const double y, const double z) const { double gLevelsetPopcorn::operator() (const double x, const double y, const double z) const
{
double s2 = (sigma)*(sigma); double s2 = (sigma)*(sigma);
double r = sqrt((x-xc)*(x-xc)+(y-yc)*(y-yc)+(z-zc)*(z-zc)); double r = sqrt((x-xc)*(x-xc)+(y-yc)*(y-yc)+(z-zc)*(z-zc));
//printf("z=%g zc=%g r=%g \n", z, zc, r); //printf("z=%g zc=%g r=%g \n", z, zc, r);
...@@ -718,7 +768,9 @@ double gLevelsetPopcorn::operator() (const double x, const double y, const doubl ...@@ -718,7 +768,9 @@ double gLevelsetPopcorn::operator() (const double x, const double y, const doubl
return val; return val;
} }
gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag) : gLevelsetPrimitive(tag) { gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag)
: gLevelsetPrimitive(tag)
{
std::vector<std::string> expressions(1, f); std::vector<std::string> expressions(1, f);
std::vector<std::string> variables(3); std::vector<std::string> variables(3);
variables[0] = "x"; variables[0] = "x";
...@@ -727,7 +779,8 @@ gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag) : gLevelsetPrimitiv ...@@ -727,7 +779,8 @@ gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag) : gLevelsetPrimitiv
_expr = new mathEvaluator(expressions, variables); _expr = new mathEvaluator(expressions, variables);
} }
double gLevelsetMathEval::operator() (const double x, const double y, const double z) const { double gLevelsetMathEval::operator() (const double x, const double y, const double z) const
{
std::vector<double> values(3), res(1); std::vector<double> values(3), res(1);
values[0] = x; values[0] = x;
values[1] = y; values[1] = y;
...@@ -736,8 +789,9 @@ double gLevelsetMathEval::operator() (const double x, const double y, const doub ...@@ -736,8 +789,9 @@ double gLevelsetMathEval::operator() (const double x, const double y, const doub
return 1.; return 1.;
} }
gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag) : gLevelsetPrimitive(tag), _box(NULL) { gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
: gLevelsetPrimitive(tag), _box(NULL)
{
modelBox = new GModel(); modelBox = new GModel();
modelBox->load(box+std::string(".msh")); modelBox->load(box+std::string(".msh"));
modelGeom = new GModel(); modelGeom = new GModel();
...@@ -745,7 +799,6 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag) ...@@ -745,7 +799,6 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
//EMI FIXME THIS //EMI FIXME THIS
int levels = 3; int levels = 3;
int refineCurvedSurfaces = 0;
// double rmax = 0.1; // double rmax = 0.1;
// double sampling = std::min(rmax, std::min(lx, std::min(ly, lz))); // double sampling = std::min(rmax, std::min(lx, std::min(ly, lz)));
...@@ -846,17 +899,18 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag) ...@@ -846,17 +899,18 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
_box->renumberNodes(); _box->renumberNodes();
_box->writeMSH("yeah.msh", false); _box->writeMSH("yeah.msh", false);
} }
double gLevelsetDistGeom::operator() (const double x, const double y, const double z) const { double gLevelsetDistGeom::operator() (const double x, const double y, const double z) const
{
double dist = _box->getValueContainingPoint(x,y,z); double dist = _box->getValueContainingPoint(x,y,z);
return dist; return dist;
} }
#if defined (HAVE_POST) #if defined (HAVE_POST)
gLevelsetPostView::gLevelsetPostView(int index, int tag) : gLevelsetPrimitive(tag), _viewIndex(index){ gLevelsetPostView::gLevelsetPostView(int index, int tag)
: gLevelsetPrimitive(tag), _viewIndex(index)
{
if(_viewIndex >= 0 && _viewIndex < (int)PView::list.size()){ if(_viewIndex >= 0 && _viewIndex < (int)PView::list.size()){
PView *view = PView::list[_viewIndex]; PView *view = PView::list[_viewIndex];
_octree = new OctreePost(view); _octree = new OctreePost(view);
...@@ -867,7 +921,8 @@ gLevelsetPostView::gLevelsetPostView(int index, int tag) : gLevelsetPrimitive(ta ...@@ -867,7 +921,8 @@ gLevelsetPostView::gLevelsetPostView(int index, int tag) : gLevelsetPrimitive(ta
} }
} }
double gLevelsetPostView::operator () (const double x, const double y, const double z) const { double gLevelsetPostView::operator () (const double x, const double y, const double z) const
{
if(!_octree) return 1.; if(!_octree) return 1.;
double val = 1.; double val = 1.;
_octree->searchScalar(x, y, z, &val, 0); _octree->searchScalar(x, y, z, &val, 0);
...@@ -875,8 +930,10 @@ double gLevelsetPostView::operator () (const double x, const double y, const dou ...@@ -875,8 +930,10 @@ double gLevelsetPostView::operator () (const double x, const double y, const dou
} }
#endif #endif
gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir, const double &R, gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir,
int tag) : gLevelsetQuadric(tag) { const double &R, int tag)
: gLevelsetQuadric(tag)
{
A[0][0] = 1.; A[0][0] = 1.;
A[1][1] = 1.; A[1][1] = 1.;
C = - R * R; C = - R * R;
...@@ -885,10 +942,14 @@ gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir, ...@@ -885,10 +942,14 @@ gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir,
rotate(rot); rotate(rot);
translate(pt); translate(pt);
} }
gLevelsetGenCylinder::gLevelsetGenCylinder (const gLevelsetGenCylinder& lv) : gLevelsetQuadric(lv){}
gLevelsetGenCylinder::gLevelsetGenCylinder (const gLevelsetGenCylinder& lv)
: gLevelsetQuadric(lv){}
gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, const double &a, gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, const double &a,
const double &b, const double &c, int tag) : gLevelsetQuadric(tag) { const double &b, const double &c, int tag)
: gLevelsetQuadric(tag)
{
A[0][0] = 1. / (a * a); A[0][0] = 1. / (a * a);
A[1][1] = 1. / (b * b); A[1][1] = 1. / (b * b);
A[2][2] = 1. / (c * c); A[2][2] = 1. / (c * c);
...@@ -898,9 +959,13 @@ gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, cons ...@@ -898,9 +959,13 @@ gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, cons
rotate(rot); rotate(rot);
translate(pt); translate(pt);
} }
gLevelsetEllipsoid::gLevelsetEllipsoid (const gLevelsetEllipsoid& lv) : gLevelsetQuadric(lv){}
gLevelsetCone::gLevelsetCone(const double *pt, const double *dir, const double &angle, int tag) : gLevelsetQuadric(tag) { gLevelsetEllipsoid::gLevelsetEllipsoid (const gLevelsetEllipsoid& lv)
: gLevelsetQuadric(lv){}
gLevelsetCone::gLevelsetCone(const double *pt, const double *dir, const double &angle,
int tag) : gLevelsetQuadric(tag)
{
A[0][0] = 1.; A[0][0] = 1.;
A[1][1] = 1.; A[1][1] = 1.;
A[2][2] = -tan(angle) * tan(angle); A[2][2] = -tan(angle) * tan(angle);
...@@ -910,10 +975,15 @@ gLevelsetCone::gLevelsetCone(const double *pt, const double *dir, const double & ...@@ -910,10 +975,15 @@ gLevelsetCone::gLevelsetCone(const double *pt, const double *dir, const double &
translate(pt); translate(pt);
} }
gLevelsetCone::gLevelsetCone (const gLevelsetCone& lv) : gLevelsetQuadric(lv) gLevelsetCone::gLevelsetCone (const gLevelsetCone& lv)
{} : gLevelsetQuadric(lv){}
gLevelsetGeneralQuadric::gLevelsetGeneralQuadric(const double *pt, const double *dir, const double &x2, const double &y2, const double &z2,
const double &z, const double &c, int tag) : gLevelsetQuadric(tag) { gLevelsetGeneralQuadric::gLevelsetGeneralQuadric(const double *pt, const double *dir,
const double &x2, const double &y2,
const double &z2, const double &z,
const double &c, int tag)
: gLevelsetQuadric(tag)
{
A[0][0] = x2; A[0][0] = x2;
A[1][1] = y2; A[1][1] = y2;
A[2][2] = z2; A[2][2] = z2;
...@@ -925,8 +995,8 @@ gLevelsetGeneralQuadric::gLevelsetGeneralQuadric(const double *pt, const double ...@@ -925,8 +995,8 @@ gLevelsetGeneralQuadric::gLevelsetGeneralQuadric(const double *pt, const double
translate(pt); translate(pt);
} }
gLevelsetGeneralQuadric::gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& lv) : gLevelsetQuadric(lv) gLevelsetGeneralQuadric::gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& lv)
{} : gLevelsetQuadric(lv){}
gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv) : gLevelset(lv) gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv) : gLevelset(lv)
{ {
...@@ -936,19 +1006,26 @@ gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv) : gLevelset(lv) ...@@ -936,19 +1006,26 @@ gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv) : gLevelset(lv)
for(unsigned i = 0; i < siz; ++i) for(unsigned i = 0; i < siz; ++i)
children[i] = _children[i]->clone(); children[i] = _children[i]->clone();
} }
gLevelsetImproved::gLevelsetImproved(const gLevelsetImproved &lv) : gLevelset(lv){
gLevelsetImproved::gLevelsetImproved(const gLevelsetImproved &lv)
: gLevelset(lv)
{
Ls = lv.Ls->clone(); Ls = lv.Ls->clone();
} }
gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *dir2, const double *dir3, gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *dir2,
const double &a, const double &b, const double &c, int tag) : gLevelsetImproved() { const double *dir3, const double &a, const double &b,
const double &c, int tag)
: gLevelsetImproved()
{
double dir1m[3] = {-dir1[0], -dir1[1], -dir1[2]}; double dir1m[3] = {-dir1[0], -dir1[1], -dir1[2]};
double dir2m[3] = {-dir2[0], -dir2[1], -dir2[2]}; double dir2m[3] = {-dir2[0], -dir2[1], -dir2[2]};
double dir3m[3] = {-dir3[0], -dir3[1], -dir3[2]}; double dir3m[3] = {-dir3[0], -dir3[1], -dir3[2]};
double n1[3]; norm(dir1, n1); double n1[3]; norm(dir1, n1);
double n2[3]; norm(dir2, n2); double n2[3]; norm(dir2, n2);
double n3[3]; norm(dir3, n3); double n3[3]; norm(dir3, n3);
double pt2[3] = {pt[0] + a * n1[0] + b * n2[0] + c * n3[0], pt[1] + a * n1[1] + b * n2[1] + c * n3[1], double pt2[3] = {pt[0] + a * n1[0] + b * n2[0] + c * n3[0],
pt[1] + a * n1[1] + b * n2[1] + c * n3[1],
pt[2] + a * n1[2] + b * n2[2] + c * n3[2]}; pt[2] + a * n1[2] + b * n2[2] + c * n3[2]};
std::vector<gLevelset *> p; std::vector<gLevelset *> p;
p.push_back(new gLevelsetPlane(pt2, dir3, tag++)); p.push_back(new gLevelsetPlane(pt2, dir3, tag++));
...@@ -960,13 +1037,18 @@ gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *d ...@@ -960,13 +1037,18 @@ gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *d
Ls = new gLevelsetIntersection(p); Ls = new gLevelsetIntersection(p);
} }
gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *pt3, const double *pt4, gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *pt3,
const double *pt5, const double *pt6, const double *pt7, const double *pt8, int tag) : gLevelsetImproved() { const double *pt4, const double *pt5, const double *pt6,
if(!isPlanar(pt1, pt2, pt3, pt4) || !isPlanar(pt5, pt6, pt7, pt8) || !isPlanar(pt1, pt2, pt5, pt6) || const double *pt7, const double *pt8, int tag)
!isPlanar(pt3, pt4, pt7, pt8) || !isPlanar(pt1, pt4, pt5, pt8) || !isPlanar(pt2, pt3, pt6, pt7)) : gLevelsetImproved()
{
if(!isPlanar(pt1, pt2, pt3, pt4) || !isPlanar(pt5, pt6, pt7, pt8) ||
!isPlanar(pt1, pt2, pt5, pt6) || !isPlanar(pt3, pt4, pt7, pt8) ||
!isPlanar(pt1, pt4, pt5, pt8) || !isPlanar(pt2, pt3, pt6, pt7))
printf("WARNING : faces of the box are not planar! %d, %d, %d, %d, %d, %d\n", printf("WARNING : faces of the box are not planar! %d, %d, %d, %d, %d, %d\n",
isPlanar(pt1, pt2, pt3, pt4), isPlanar(pt5, pt6, pt7, pt8), isPlanar(pt1, pt2, pt5, pt6), isPlanar(pt1, pt2, pt3, pt4), isPlanar(pt5, pt6, pt7, pt8),
isPlanar(pt3, pt4, pt7, pt8), isPlanar(pt1, pt4, pt5, pt8), isPlanar(pt2, pt3, pt6, pt7)); isPlanar(pt1, pt2, pt5, pt6), isPlanar(pt3, pt4, pt7, pt8),
isPlanar(pt1, pt4, pt5, pt8), isPlanar(pt2, pt3, pt6, pt7));
std::vector<gLevelset *> p; std::vector<gLevelset *> p;
p.push_back(new gLevelsetPlane(pt5, pt6, pt8, tag++)); p.push_back(new gLevelsetPlane(pt5, pt6, pt8, tag++));
p.push_back(new gLevelsetPlane(pt1, pt4, pt2, tag++)); p.push_back(new gLevelsetPlane(pt1, pt4, pt2, tag++));
...@@ -977,9 +1059,14 @@ gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *p ...@@ -977,9 +1059,14 @@ gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *p
Ls = new gLevelsetIntersection(p); Ls = new gLevelsetIntersection(p);
} }
gLevelsetBox::gLevelsetBox(const gLevelsetBox &lv) : gLevelsetImproved(lv){} gLevelsetBox::gLevelsetBox(const gLevelsetBox &lv)
: gLevelsetImproved(lv){}
gLevelsetCylinder::gLevelsetCylinder(const std::vector<double> &pt, const std::vector<double> &dir, const double &R, const double &H, int tag) : gLevelsetImproved() { gLevelsetCylinder::gLevelsetCylinder(const std::vector<double> &pt,
const std::vector<double> &dir,
const double &R, const double &H, int tag)
: gLevelsetImproved()
{
double pt1[3]={pt[0], pt[1], pt[2]}; double pt1[3]={pt[0], pt[1], pt[2]};
double dir1[3] = {dir[0], dir[1], dir[2]}; double dir1[3] = {dir[0], dir[1], dir[2]};
double dir2[3] = {-dir1[0], -dir1[1], -dir1[2]}; double dir2[3] = {-dir1[0], -dir1[1], -dir1[2]};
...@@ -992,7 +1079,10 @@ gLevelsetCylinder::gLevelsetCylinder(const std::vector<double> &pt, const std::v ...@@ -992,7 +1079,10 @@ gLevelsetCylinder::gLevelsetCylinder(const std::vector<double> &pt, const std::v
Ls = new gLevelsetIntersection(p); Ls = new gLevelsetIntersection(p);
} }
gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir, const double &R, const double &H, int tag) : gLevelsetImproved() { gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir,
const double &R, const double &H, int tag)
: gLevelsetImproved()
{
double dir2[3] = {-dir[0], -dir[1], -dir[2]}; double dir2[3] = {-dir[0], -dir[1], -dir[2]};
double n[3]; norm(dir, n); double n[3]; norm(dir, n);
double pt2[3] = {pt[0] + H * n[0], pt[1] + H * n[1], pt[2] + H * n[2]}; double pt2[3] = {pt[0] + H * n[0], pt[1] + H * n[1], pt[2] + H * n[2]};
...@@ -1003,7 +1093,11 @@ gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir, const ...@@ -1003,7 +1093,11 @@ gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir, const
Ls = new gLevelsetIntersection(p); Ls = new gLevelsetIntersection(p);
} }
gLevelsetCylinder::gLevelsetCylinder(const double * pt, const double *dir, const double &R, const double &r, const double &H, int tag) : gLevelsetImproved() { gLevelsetCylinder::gLevelsetCylinder(const double * pt, const double *dir,
const double &R, const double &r, const double &H,
int tag)
: gLevelsetImproved()
{
double dir2[3] = {-dir[0], -dir[1], -dir[2]}; double dir2[3] = {-dir[0], -dir[1], -dir[2]};
double n[3]; norm(dir, n); double n[3]; norm(dir, n);
double pt2[3] = {pt[0] + H * n[0], pt[1] + H * n[1], pt[2] + H * n[2]}; double pt2[3] = {pt[0] + H * n[0], pt[1] + H * n[1], pt[2] + H * n[2]};
...@@ -1016,20 +1110,30 @@ gLevelsetCylinder::gLevelsetCylinder(const double * pt, const double *dir, const ...@@ -1016,20 +1110,30 @@ gLevelsetCylinder::gLevelsetCylinder(const double * pt, const double *dir, const
p2.push_back(new gLevelsetGenCylinder(pt, dir, r, tag)); p2.push_back(new gLevelsetGenCylinder(pt, dir, r, tag));
Ls = new gLevelsetCut(p2); Ls = new gLevelsetCut(p2);
} }
gLevelsetCylinder::gLevelsetCylinder(const gLevelsetCylinder &lv) : gLevelsetImproved(lv){}
gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1, const double *dir2, gLevelsetCylinder::gLevelsetCylinder(const gLevelsetCylinder &lv)
const double &H1, const double &H2, const double &H3, : gLevelsetImproved(lv){}
const double &R1, const double &r1, const double &R2, const double &r2,
const double &L1, const double &L2, const double &E, int tag) : gLevelsetImproved() { gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1,
const double *dir2, const double &H1,
const double &H2, const double &H3,
const double &R1, const double &r1,
const double &R2, const double &r2,
const double &L1, const double &L2,
const double &E, int tag)
: gLevelsetImproved()
{
double n1[3]; norm(dir1, n1); double n1[3]; norm(dir1, n1);
double n2[3]; norm(dir2, n2); double n2[3]; norm(dir2, n2);
double pt1[3] = {pt[0] - n2[0] * H1 / 2., pt[1] - n2[1] * H1 / 2., pt[2] - n2[2] * H1 / 2.}; double pt1[3] = {pt[0] - n2[0] * H1 / 2., pt[1] - n2[1] * H1 / 2.,
double pt2[3] = {pt[0] + n1[0] * E - n2[0] * H2 / 2., pt[1] + n1[1] * E - n2[1] * H2 / 2., pt[2] - n2[2] * H1 / 2.};
double pt2[3] = {pt[0] + n1[0] * E - n2[0] * H2 / 2.,
pt[1] + n1[1] * E - n2[1] * H2 / 2.,
pt[2] + n1[2] * E - n2[2] * H2 / 2.}; pt[2] + n1[2] * E - n2[2] * H2 / 2.};
double dir3[3]; cross(pt1, pt2, pt, dir3); double dir3[3]; cross(pt1, pt2, pt, dir3);
double n3[3]; norm(dir3, n3); double n3[3]; norm(dir3, n3);
double pt31[3] = {pt[0] - n2[0] * H3 / 2. + n3[0] * L1 / 2., pt[1] - n2[1] * H3 / 2. + n3[1] * L1 / 2., double pt31[3] = {pt[0] - n2[0] * H3 / 2. + n3[0] * L1 / 2.,
pt[1] - n2[1] * H3 / 2. + n3[1] * L1 / 2.,
pt[2] - n2[2] * H3 / 2. + n3[2] * L1 / 2.}; pt[2] - n2[2] * H3 / 2. + n3[2] * L1 / 2.};
double pt32[3] = {pt31[0] - n3[0] * L1, pt31[1] - n3[1] * L1, pt31[2] - n3[2] * L1}; double pt32[3] = {pt31[0] - n3[0] * L1, pt31[1] - n3[1] * L1, pt31[2] - n3[2] * L1};
double pt33[3] = {pt32[0] + n2[0] * H3, pt32[1] + n2[1] * H3, pt32[2] + n2[2] * H3}; double pt33[3] = {pt32[0] + n2[0] * H3, pt32[1] + n2[1] * H3, pt32[2] + n2[2] * H3};
...@@ -1051,4 +1155,5 @@ gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1, const dou ...@@ -1051,4 +1155,5 @@ gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1, const dou
Ls = new gLevelsetCut(p2); Ls = new gLevelsetCut(p2);
} }
gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv) : gLevelsetImproved(lv){} gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv)
: gLevelsetImproved(lv){}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment