diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index 29999e3df33ae3790d1faf7e6dbfea71ffdeec25..19391052b31de6a812db57bda1d5993011e38207 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -61,7 +61,7 @@ static void printNodes(std::set<MVertex *> myNodes)
{
FILE * xyz = fopen("myNodes.pos","w");
fprintf(xyz,"View \"\"{\n");
- 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 *v = *itv;
fprintf(xyz,"SP(%g,%g,%g){%d};\n", v->x(), v->y(), v->z(), v->getNum());
}
@@ -86,8 +86,9 @@ static void exportParametrizedMesh(fullMatrix<double> &UV, int nbNodes)
GRbf::GRbf(double sizeBox, int variableEps, int rbfFun,
std::map<MVertex*, SVector3> _normals,
std::set<MVertex *> allNodes, std::vector<MVertex*> bcNodes, bool _isLocal)
- : sBox(sizeBox), variableShapeParam(variableEps), radialFunctionIndex (rbfFun),
- _inUV(0), isLocal(_isLocal)
+ : isLocal(_isLocal), _inUV(0), sBox(sizeBox), variableShapeParam(variableEps),
+ radialFunctionIndex (rbfFun)
+
{
#if defined (HAVE_ANN)
XYZkdtree = 0;
@@ -125,7 +126,6 @@ GRbf::GRbf(double sizeBox, int variableEps, int rbfFun,
normals.resize(nbNodes,3);
int index = 0;
double dist_min = 1.e6;
- double dist_max = 1.e-6;
for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end(); ++itv){
MVertex *v1 = *itv;
centers(index,0) = v1->x()/sBox;
@@ -143,7 +143,7 @@ GRbf::GRbf(double sizeBox, int variableEps, int rbfFun,
double dist = sqrt((v1->x()-v2->x())*(v1->x()-v2->x())+
(v1->y()-v2->y())*(v1->y()-v2->y())+
(v1->z()-v2->z())*(v1->z()-v2->z()))/sBox;
- if (dist<dist_min && *itv != *itv2 && dist > 1.e-5) dist_min = dist;
+ if (dist < dist_min && *itv != *itv2 && dist > 1.e-5) dist_min = dist;
}
index++;
}
@@ -196,7 +196,7 @@ void GRbf::buildOctree(double radius)
{
//printf("building octree radius = %g \n", radius);
SBoundingBox3d bb;
- for (int i= 0; i< nbNodes; i++)
+ for (int i = 0; i < nbNodes; i++)
bb += SPoint3(centers(i,0),centers(i,1), centers(i,2));
double origin[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
double ssize[3] = {bb.max().x() - bb.min().x(),
@@ -207,7 +207,7 @@ void GRbf::buildOctree(double radius)
SphereCentroid, SphereInEle);
Sphere *_sph = new Sphere[nbNodes];
- for (int i= 0; i< nbNodes; i++){
+ for (int i = 0; i < nbNodes; i++){
_sph[i].index = i;
_sph[i].radius = radius;
_sph[i].center = SPoint3(centers(i,0),centers(i,1), centers(i,2));
@@ -215,7 +215,7 @@ void GRbf::buildOctree(double radius)
}
Octree_Arrange(oct);
- for (int i= 0; i< nbNodes; i++){
+ for (int i = 0; i < nbNodes; i++){
std::vector<void*> l;
double P[3] = {centers(i,0),centers(i,1), centers(i,2)};
Octree_SearchAll(P, oct, &l);
@@ -223,8 +223,7 @@ void GRbf::buildOctree(double radius)
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++) {
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());
}
@@ -279,13 +278,13 @@ void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
if(nodesInSphere.size() == 0) buildOctree(radius);
fullMatrix<double> curvature(cntrs.size1(), 1);
- for (int i=0;i<numNodes ; ++i) {
+ for (int i = 0; i < numNodes; ++i) {
std::vector<int> &pts = nodesInSphere[i];
fullMatrix<double> nodes_in_sph(pts.size(),3);
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,1) = cntrs(pts[k],1);
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)
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,
@@ -416,7 +416,7 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
fullMatrix<double> sz(numNodes,1),norms(numNodes,3), cntrsPlus(numNodes+1,3);
//Computes the normal vectors to the surface at each node
- for (int i=0;i<numNodes ; ++i){
+ for (int i = 0; i < numNodes; ++i){
ONES(i,0) = 0.0;
cntrsPlus(i,0) = cntrs(i,0);
cntrsPlus(i,1) = cntrs(i,1);
@@ -430,7 +430,7 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
evalRbfDer(1,cntrsPlus,cntrs,ONES,sx);
evalRbfDer(2,cntrsPlus,cntrs,ONES,sy);
evalRbfDer(3,cntrsPlus,cntrs,ONES,sz);
- for (int i=0;i<numNodes ; ++i){
+ for (int i = 0; i < numNodes; ++i){
//GMSH NORMALS
//sx(i,0) = normals(i,0);
//sy(i,0) = normals(i,1);
@@ -443,8 +443,8 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
norms(i,0) = sx(i,0);norms(i,1) = sy(i,0);norms(i,2) = sz(i,0);
}
- for (int i=0;i<numNodes ; ++i){
- for (int j=0;j<3 ; ++j){
+ for (int i = 0; i < numNodes; ++i){
+ for (int j = 0; j < 3; ++j){
level_set_nodes(i,j) = cntrs(i,j);
level_set_nodes(i+numNodes,j) = cntrs(i,j)-delta*norms(i,j);
level_set_nodes(i+2*numNodes,j) = cntrs(i,j)+delta*norms(i,j);
@@ -460,8 +460,8 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
}
void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs,
- const fullMatrix<double> &normals,
- fullMatrix<double> &Oper)
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &Oper)
{
isLocal = true;
int numNodes = cntrs.size1();
@@ -469,7 +469,7 @@ void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs,
if(nodesInSphere.size() == 0) buildOctree(radius);
- for (int i=0;i<numNodes ; ++i){
+ for (int i = 0; i < numNodes; ++i){
std::vector<int> &pts = nodesInSphere[i];
fullMatrix<double> nodes_in_sph(pts.size(),3), local_normals(pts.size(),3);
@@ -477,23 +477,23 @@ void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs,
LocalOper.setAll(0.0);
- for (int k=0; k< pts.size(); k++){
- nodes_in_sph(k,0) = cntrs(pts[k],0);
- nodes_in_sph(k,1) = cntrs(pts[k],1);
- nodes_in_sph(k,2) = cntrs(pts[k],2);
- local_normals(k,0)=normals(pts[k],0);
- local_normals(k,1)=normals(pts[k],1);
- local_normals(k,2)=normals(pts[k],2);
+ for (unsigned int k = 0; k < pts.size(); k++){
+ nodes_in_sph(k, 0) = cntrs(pts[k], 0);
+ nodes_in_sph(k, 1) = cntrs(pts[k], 1);
+ nodes_in_sph(k, 2) = cntrs(pts[k], 2);
+ local_normals(k, 0) = normals(pts[k], 0);
+ local_normals(k, 1) = normals(pts[k], 1);
+ local_normals(k, 2) = normals(pts[k], 2);
}
RbfLapSurface_global_projection(nodes_in_sph,local_normals, LocalOper);
- for (int j=0;j<pts.size() ; j++)
- Oper(i,pts[j])=LocalOper(0,j);
+ for (unsigned int j = 0; j < pts.size(); j++)
+ Oper(i, pts[j]) = LocalOper(0, j);
}
}
-void GRbf::RbfLapSurface_global_projection( const fullMatrix<double> &cntrs,
+void GRbf::RbfLapSurface_global_projection(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
fullMatrix<double> &Oper)
{
@@ -513,7 +513,7 @@ void GRbf::RbfLapSurface_global_projection( const fullMatrix<double> &cntrs,
// Normalizes
double norm;
- for (int i=0;i<numNodes;i++){
+ for (int i = 0; i < numNodes;i++){
norm = sqrt(sx(i,0)*sx(i,0)+sy(i,0)*sy(i,0)+sz(i,0)*sz(i,0));
sx(i,0) /= norm;
sy(i,0) /= norm;
@@ -526,8 +526,8 @@ void GRbf::RbfLapSurface_global_projection( const fullMatrix<double> &cntrs,
RbfOp(3,cntrs,cntrs,Dz);
// Fills up the operator matrix
- for (int i=0;i<numNodes ; ++i){
- for (int j=0;j<numNodes ; ++j){
+ for (int i = 0; i < numNodes; ++i){
+ for (int j = 0; j < numNodes; ++j){
PDx(i,j) = (1-sx(i,0)*sx(i,0))*Dx(i,j)-sx(i,0)*sy(i,0)*Dy(i,j)-sx(i,0)*sz(i,0)*Dz(i,j);
PDy(i,j) = -sx(i,0)*sy(i,0)*Dx(i,j)+(1-sy(i,0)*sy(i,0))*Dy(i,j)-sy(i,0)*sz(i,0)*Dz(i,j);
PDz(i,j) = -sx(i,0)*sz(i,0)*Dx(i,j)-sy(i,0)*sz(i,0)*Dy(i,j)+(1-sz(i,0)*sz(i,0))*Dz(i,j);
@@ -536,17 +536,17 @@ void GRbf::RbfLapSurface_global_projection( const fullMatrix<double> &cntrs,
PDx.mult(PDx,PDxx);
PDy.mult(PDy,PDyy);
PDz.mult(PDz,PDzz);
- for (int i=0;i<numNodes ; ++i){
- for (int j=0;j<numNodes ; ++j){
+ for (int i = 0; i < numNodes; ++i){
+ for (int j = 0; j < numNodes; ++j){
Oper(i,j) = PDxx(i,j)+PDyy(i,j)+PDzz(i,j);
}
}
}
void GRbf::RbfLapSurface_local_CPM(bool isLow,
- const fullMatrix<double> &cntrs,
- const fullMatrix<double> &normals,
- fullMatrix<double> &Oper)
+ const fullMatrix<double> &cntrs,
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &Oper)
{
isLocal = true;
int numNodes = cntrs.size1();
@@ -555,7 +555,7 @@ void GRbf::RbfLapSurface_local_CPM(bool isLow,
buildOctree(radius);
setup_level_set(cntrs,normals,extendedX,surfInterp);
- for (int i=0;i<numNodes ; ++i){
+ for (int i = 0; i < numNodes; ++i){
std::vector<int> &pts = nodesInSphere[i];
int nbp = pts.size();
fullMatrix<double> nodes_in_sph(nbp,3), local_normals(nbp,3);
@@ -574,7 +574,7 @@ void GRbf::RbfLapSurface_local_CPM(bool isLow,
if (isLow) RbfLapSurface_global_CPM_low(nodes_in_sph,local_normals,LocalOper);
else RbfLapSurface_global_CPM_high_2(nodes_in_sph,local_normals,LocalOper);
- for (int j = 0; j < nbp ; j++){
+ for (int j = 0; j < nbp; j++){
Oper(i,pts[j])=LocalOper(0,j);
Oper(i,pts[j]+numNodes)=LocalOper(0,j+nbp);
Oper(i,pts[j]+2*numNodes)=LocalOper(0,j+2*nbp);
@@ -610,30 +610,30 @@ void GRbf::RbfLapSurface_local_CPM_sparse(std::vector<MVertex*> &bndVertices, bo
buildOctree(radius);
setup_level_set(cntrs,normals,extendedX,surfInterp);
- for (int i = 0 ; i < numNodes ; ++i) {
+ for (int i = 0; i < numNodes; ++i) {
std::vector<int> &pts = nodesInSphere[i];
if (bndIndices.count(i) > 0) {
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 + 2 * numNodes, pts[j]);
}
}
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 + numNodes, pts[j]);
sys.insertInSparsityPattern(i + 2 * numNodes, pts[j]);
}
}
}
- for (int i=0;i<numNodes ; ++i){
+ for (int i = 0; i < numNodes; ++i){
std::vector<int> &pts = nodesInSphere[i];
int nbp = pts.size();
fullMatrix<double> nodes_in_sph(nbp,3), local_normals(nbp,3);
fullMatrix<double> LocalOper;
- for (int k=0; k< nbp; k++){
+ for (int k = 0; k < nbp; k++){
nodes_in_sph(k,0) = cntrs(pts[k],0);
nodes_in_sph(k,1) = cntrs(pts[k],1);
nodes_in_sph(k,2) = cntrs(pts[k],2);
@@ -650,7 +650,7 @@ void GRbf::RbfLapSurface_local_CPM_sparse(std::vector<MVertex*> &bndVertices, bo
if (isBnd) {
sys.addToMatrix(i, i, 1.);
}
- for (int j = 0; j < nbp ; j++){
+ for (int j = 0; j < nbp; j++){
if (!isBnd) {
sys.addToMatrix(i, pts[j], LocalOper(0,j));
sys.addToMatrix(i, pts[j] + numNodes, LocalOper(0,j + nbp));
@@ -712,8 +712,8 @@ void GRbf::RbfLapSurface_global_CPM_high_2(const fullMatrix<double> &cntrs,
// Fills up the operator matrix
AOper.resize(nnTot, nnTot);
- for (int i=0;i<numNodes ; ++i){
- for (int j=0;j<nnTot ; ++j){
+ for (int i = 0; i < numNodes; ++i){
+ for (int j = 0; j < nnTot; ++j){
AOper(i,j) = Alap(i,j);
AOper(i+numNodes,j)=sx(i,0)*Ax(i,j)+sy(i,0)*Ay(i,j)+sz(i,0)*Az(i,j);
AOper(i+2*numNodes,j)=sx(i,0)*sx(i,0)*Axx(i,j)+sy(i,0)*sy(i,0)*Ayy(i,j)+sz(i,0)*sz(i,0)*Azz(i,j)+2*sx(i,0)*sy(i,0)*Axy(i,j)+2*sx(i,0)*sz(i,0)*Axz(i,j)+2*sy(i,0)*sz(i,0)*Ayz(i,j)+(sx(i,0)*sxx(i,0)+sy(i,0)*sxy(i,0)+sz(i,0)*sxz(i,0))*Ax(i,j)+(sx(i,0)*sxy(i,0)+sy(i,0)*syy(i,0)+sz(i,0)*syz(i,0))*Ay(i,j)+(sx(i,0)*sxz(i,0)+sy(i,0)*syz(i,0)+sz(i,0)*szz(i,0))*Az(i,j);
@@ -769,8 +769,8 @@ void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
Alap=generateRbfMat(222,extX,cntrs);
// Fills up the operator matrix
- for (int i=0;i<numNodes ; ++i){
- for (int j=0;j<nnTot ; ++j){
+ for (int i = 0; i < numNodes; ++i){
+ for (int j = 0; j < nnTot; ++j){
AOper(i,j) = A(i,j);
AOper(i+numNodes,j)=sx(i,0)*Ax(i,j)+sy(i,0)*Ay(i,j)+sz(i,0)*Az(i,j);
AOper(i+2*numNodes,j)=sx(i,0)*sx(i,0)*Axx(i,j)+sy(i,0)*sy(i,0)*Ayy(i,j)+sz(i,0)*sz(i,0)*Azz(i,j)+2*sx(i,0)*sy(i,0)*Axy(i,j)+2*sx(i,0)*sz(i,0)*Axz(i,j)+2*sy(i,0)*sz(i,0)*Ayz(i,j)+(sx(i,0)*sxx(i,0)+sy(i,0)*sxy(i,0)+sz(i,0)*sxz(i,0))*Ax(i,j)+(sx(i,0)*sxy(i,0)+sy(i,0)*syy(i,0)+sz(i,0)*syz(i,0))*Ay(i,j)+(sx(i,0)*sxz(i,0)+sy(i,0)*syz(i,0)+sz(i,0)*szz(i,0))*Az(i,j);
@@ -778,8 +778,8 @@ void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
}
AOper.invertInPlace();
Alap.mult(AOper,Temp);
- for (int i=0;i<numNodes ; ++i){
- for (int j=0;j<numNodes ; ++j){
+ for (int i = 0; i < numNodes; ++i){
+ for (int j = 0; j < numNodes; ++j){
Oper(i,j) = Temp(i,j);
}
}
@@ -788,8 +788,8 @@ void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
}
void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs,
- const fullMatrix<double> &normals,
- fullMatrix<double> &Oper)
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &Oper)
{
int numNodes = cntrs.size1();
int nnTot = 3*numNodes;
@@ -808,7 +808,7 @@ void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs,
evalRbfDer(1,extX,extX,surf,sx);
evalRbfDer(2,extX,extX,surf,sy);
evalRbfDer(3,extX,extX,surf,sz);
- for (int i=0;i<nnTot ; ++i){
+ for (int i = 0; i < nnTot; ++i){
double normFactor = sqrt(sx(i,0)*sx(i,0)+sy(i,0)*sy(i,0)+sz(i,0)*sz(i,0));
sx(i,0) = sx(i,0)/normFactor;
sy(i,0) = sy(i,0)/normFactor;
@@ -817,8 +817,8 @@ void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs,
}
//Computes the inside and outside layers nodes of the surface
- for (int i=0;i<numNodes ; ++i){
- for (int j=0;j<3 ; ++j){
+ for (int i = 0; i < numNodes; ++i){
+ for (int j = 0; j < 3; ++j){
extRBFX(i,j) = cntrs(i,j)-delta*norm(i,j);
extRBFX(i+numNodes,j) = cntrs(i,j)+delta*norm(i,j);
}
@@ -829,12 +829,12 @@ void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs,
RbfOp(222,extX,cntrs,PLap); //'222' for the Laplacian
// Fills up the operator matrix
- for (int i=0; i < numNodes; i++){
- for (int j=0; j < nnTot ; ++j){
+ for (int i = 0; i < numNodes; i++){
+ for (int j = 0; j < nnTot; ++j){
Oper(i,j) = PLap(i,j);
double del = (i == 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 pos1 = (i+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+2*numNodes,j) = del + Ix2extX(i+numNodes,j); //+ pos2; //
}
@@ -843,9 +843,9 @@ void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs,
}
void GRbf::solveHarmonicMap_sparse(linearSystem<double> &sys, int numNodes,
- const std::vector<MVertex*> &bcNodes,
- const std::vector<double> &coords,
- std::map<MVertex*, SPoint3> &rbf_param)
+ const std::vector<MVertex*> &bcNodes,
+ const std::vector<double> &coords,
+ std::map<MVertex*, SPoint3> &rbf_param)
{
#if defined(HAVE_SOLVER)
Msg::Info("*** RBF ... solving map");
@@ -854,7 +854,7 @@ void GRbf::solveHarmonicMap_sparse(linearSystem<double> &sys, int numNodes,
for (int j = 0; j < 2; ++j) {
sys.zeroRightHandSide();
//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]);
if (itN != myNodes.end()){
std::map<MVertex*, int>::iterator itm = _mapV.find(bcNodes[i]);
@@ -870,12 +870,12 @@ void GRbf::solveHarmonicMap_sparse(linearSystem<double> &sys, int numNodes,
}
//SQUARE
- // 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 *v = *itv;
// if (v->getNum() == 1){ //2900){
// std::map<MVertex*, int>::iterator itm = _mapV.find(v);
// int iFix = itm->second;
- // for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
+ // for (int j = 0; j < nb; ++j) Oper(iFix,j) = 0.0;
// Oper(iFix,iFix) = 1.0;
// F(iFix,0) = 0.0;
// F(iFix,1) = 0.0;
@@ -883,7 +883,7 @@ void GRbf::solveHarmonicMap_sparse(linearSystem<double> &sys, int numNodes,
// if (v->getNum() == 2){//1911){
// std::map<MVertex*, int>::iterator itm = _mapV.find(v);
// int iFix = itm->second;
- // for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
+ // for (int j = 0; j < nb; ++j) Oper(iFix,j) = 0.0;
// Oper(iFix,iFix) = 1.0;
// F(iFix,0) = 1.0;
// F(iFix,1) = 0.0;
@@ -891,7 +891,7 @@ void GRbf::solveHarmonicMap_sparse(linearSystem<double> &sys, int numNodes,
// if (v->getNum() == 3){//4844){
// std::map<MVertex*, int>::iterator itm = _mapV.find(v);
// int iFix = itm->second;
- // for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
+ // for (int j = 0; j < nb; ++j) Oper(iFix,j) = 0.0;
// Oper(iFix,iFix) = 1.0;
// F(iFix,0) = 1.0;
// F(iFix,1) = 1.0;
@@ -899,7 +899,7 @@ void GRbf::solveHarmonicMap_sparse(linearSystem<double> &sys, int numNodes,
// if (v->getNum() == 4){//3187){
// std::map<MVertex*, int>::iterator itm = _mapV.find(v);
// int iFix = itm->second;
- // for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
+ // for (int j = 0; j < nb; ++j) Oper(iFix,j) = 0.0;
// Oper(iFix,iFix) = 1.0;
// F(iFix,0) = 0.0;
// F(iFix,1) = 1.0;
@@ -954,13 +954,13 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
F.scale(0.0);
//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]);
if (itN != myNodes.end()){
std::map<MVertex*, int>::iterator itm = _mapV.find(bcNodes[i]);
double theta = 2 * M_PI * coords[i];
int iFix = itm->second;
- for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
+ for (int j = 0; j < nb; ++j) Oper(iFix,j) = 0.0;
Oper(iFix,iFix) = 1.0;
F(iFix,0) = cos(theta);
F(iFix,1) = sin(theta);
@@ -968,12 +968,12 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
}
//SQUARE
- // 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 *v = *itv;
// if (v->getNum() == 1){ //2900){
// std::map<MVertex*, int>::iterator itm = _mapV.find(v);
// int iFix = itm->second;
- // for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
+ // for (int j = 0; j < nb; ++j) Oper(iFix,j) = 0.0;
// Oper(iFix,iFix) = 1.0;
// F(iFix,0) = 0.0;
// F(iFix,1) = 0.0;
@@ -981,7 +981,7 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
// if (v->getNum() == 2){//1911){
// std::map<MVertex*, int>::iterator itm = _mapV.find(v);
// int iFix = itm->second;
- // for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
+ // for (int j = 0; j < nb; ++j) Oper(iFix,j) = 0.0;
// Oper(iFix,iFix) = 1.0;
// F(iFix,0) = 1.0;
// F(iFix,1) = 0.0;
@@ -989,7 +989,7 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
// if (v->getNum() == 3){//4844){
// std::map<MVertex*, int>::iterator itm = _mapV.find(v);
// int iFix = itm->second;
- // for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
+ // for (int j = 0; j < nb; ++j) Oper(iFix,j) = 0.0;
// Oper(iFix,iFix) = 1.0;
// F(iFix,0) = 1.0;
// F(iFix,1) = 1.0;
@@ -997,7 +997,7 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
// if (v->getNum() == 4){//3187){
// std::map<MVertex*, int>::iterator itm = _mapV.find(v);
// int iFix = itm->second;
- // for (int j=0;j<nb ; ++j) Oper(iFix,j) = 0.0;
+ // for (int j = 0; j < nb; ++j) Oper(iFix,j) = 0.0;
// Oper(iFix,iFix) = 1.0;
// F(iFix,0) = 0.0;
// F(iFix,1) = 1.0;
@@ -1108,163 +1108,7 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
lastX = XX;
lastY = YY;
lastZ = ZZ;
- lastDXDU = dXdu ;
- lastDXDV = dXdv ;
+ lastDXDU = dXdu;
+ lastDXDV = dXdv;
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
-// }
-
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index f0515ed706f0e1a1807b4e7cf1b76ae1a6387bbf..3cf149554a9eae9835008c788d6d66b4d82643f8 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -31,6 +31,7 @@ void insertActiveCells(double x, double y, double z, double rmax,
for (int k = k1; k <= k2; k++)
box.insertActiveCell(box.getCellIndex(i, j, k));
}
+
void fillPointCloud(GEdge *ge, double sampling, std::vector<SPoint3> &points)
{
Range<double> t_bounds = ge->parBounds(0);
@@ -71,7 +72,7 @@ int removeBadChildCells(cartesianBox<double> *parent)
if(exist[ii]) atLeastOne = true;
if(!exist[ii]) butNotAll = true;
}
- if(atLeastOne &&
+ if(atLeastOne &&
(butNotAll ||
(i != 0 && !parent->activeCellExists(parent->getCellIndex(i - 1, j, k))) ||
(i != I - 1 && !parent->activeCellExists(parent->getCellIndex(i + 1, j, k))) ||
@@ -108,11 +109,9 @@ void removeParentCellsWithChildren(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<int> indices;
- for (cartesianBox<double>::valIter it = box.nodalValuesBegin();
+ for (cartesianBox<double>::valIter it = box.nodalValuesBegin();
it != box.nodalValuesEnd(); ++it){
nodes.push_back(box.getNodeCoordinates(it->first));
indices.push_back(it->first);
@@ -122,7 +121,7 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
std::vector<SPoint3> dummy;
for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
if((*fit)->geomType() == GEntity::DiscreteSurface){
- for(unsigned int k = 0; k < (*fit)->getNumMeshElements(); k++){
+ for(unsigned int k = 0; k < (*fit)->getNumMeshElements(); k++){
std::vector<double> iDistances;
std::vector<SPoint3> iClosePts;
std::vector<double> iDistancesE;
@@ -137,9 +136,9 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
//sign plus if in the direction of the normal
signedDistancesPointsTriangle(localdist, dummy, nodes, p2, p1, p3);
}
- if(dist.empty())
+ if(dist.empty())
dist = localdist;
- else{
+ else{
for (unsigned int j = 0; j < localdist.size(); j++){
dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
}
@@ -160,21 +159,24 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
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) {
- return d11 * (d22 * d33 - d23 * d32) - d21 * (d12 * d33 - d13 * d32) + d31 * (d12 * d23 - d13 * d22);
+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);
}
-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]);
norm[0] = vec[0] / n;
norm[1] = vec[1] / 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 v2[3] = {pt2[0] - pt0[0], pt2[1] - pt0[1], pt2[2] - pt0[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
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 c2[3]; cross(pt1, pt2, pt4, c2);
double n1[3]; norm(c1, n1);
@@ -190,9 +194,10 @@ 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]);
}
-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
+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
switch (index) {
case 0 : // GA
switch (p) {
@@ -207,22 +212,26 @@ inline double evalRadialFnDer (int p, int index, double dx, double dy, double dz
case 1: return ep*ep*dx/sqrt(1+ep*ep*r2);
case 2: return ep*ep*dy/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");
fprintf(xyz,"View \"\"{\n");
- 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));
+ 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,"};\n");
fclose(xyz);
}
// 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;
Q.push(this);
while(!Q.empty()){
@@ -238,8 +247,10 @@ void gLevelset::getPrimitives(std::vector<gLevelset *> &gLsPrimitives) {
}
}
}
+
// 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 *> Sc; // levelset checked
S.push(this);
@@ -268,7 +279,8 @@ void gLevelset::getPrimitivesPO(std::vector<gLevelset *> &gLsPrimitives) {
}
// 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 *> Sc; // levelset checked
S.push(this);
@@ -302,25 +314,37 @@ gLevelset::gLevelset(const gLevelset &lv)
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];
b = norm[1];
c = norm[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];
b = norm[1];
c = norm[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]);
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.);
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;
b = lv.b;
c = lv.c;
@@ -330,12 +354,13 @@ gLevelsetPlane::gLevelsetPlane(const gLevelsetPlane &lv) : gLevelsetPrimitive(lv
//level set defined by points (RBF interpolation)
fullMatrix<double> gLevelsetPoints::generateRbfMat(int p, int index,
const fullMatrix<double> &nodes1,
- const fullMatrix<double> &nodes2) const {
+ const fullMatrix<double> &nodes2) const
+{
int m = nodes2.size1();
int n = nodes1.size1();
fullMatrix<double> rbfMat(m,n);
- double eps =0.5/delta;
+ double eps =0.5/delta;
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
double dx = nodes2(i,0)-nodes1(j,0);
@@ -347,12 +372,13 @@ fullMatrix<double> gLevelsetPoints::generateRbfMat(int p, int index,
return rbfMat;
}
-void gLevelsetPoints::RbfOp(int p, int index,
- const fullMatrix<double> &cntrs,
- const fullMatrix<double> &nodes,
- fullMatrix<double> &D,
- bool isLocal) const {
+void gLevelsetPoints::RbfOp(int p, int index,
+ const fullMatrix<double> &cntrs,
+ const fullMatrix<double> &nodes,
+ fullMatrix<double> &D,
+ bool isLocal) const
+{
fullMatrix<double> rbfMatB = generateRbfMat(p,index, cntrs,nodes);
//printf("size=%d %d \n", rbfMatB.size1(), rbfMatB.size2());
//rbfMatB.print("MatB");
@@ -373,37 +399,35 @@ void gLevelsetPoints::RbfOp(int p, int index,
}
-
void gLevelsetPoints::evalRbfDer(int p, int index,
- const fullMatrix<double> &cntrs,
- const fullMatrix<double> &nodes,
- const fullMatrix<double> &fValues,
- fullMatrix<double> &fApprox, bool isLocal) const {
-
+ const fullMatrix<double> &cntrs,
+ const fullMatrix<double> &nodes,
+ const fullMatrix<double> &fValues,
+ fullMatrix<double> &fApprox, bool isLocal) const
+{
fApprox.resize(nodes.size1(),fValues.size2());
fullMatrix<double> D;
RbfOp(p,index, cntrs,nodes,D,isLocal);
//D.print("D");
fApprox.gemm(D,fValues, 1.0, 0.0);
-
}
-
void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs,
- fullMatrix<double> &level_set_nodes,
- fullMatrix<double> &level_set_funvals){
-
+ fullMatrix<double> &level_set_nodes,
+ fullMatrix<double> &level_set_funvals)
+{
int numNodes = cntrs.size1();
int nTot = 3*numNodes;
double normFactor;
level_set_nodes.resize(nTot,3);
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
double dist_min = 1.e6;
double dist_max = 1.e-6;
- for (int i=0;i<numNodes ; ++i){
+ for (int i = 0; i < numNodes; ++i){
ONES(i,0)=1.0;
cntrsPlus(i,0) = cntrs(i,0);
cntrsPlus(i,1) = cntrs(i,1);
@@ -425,19 +449,19 @@ void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs,
evalRbfDer(1, 1, cntrsPlus,cntrs,ONES,sx, true);
evalRbfDer(2, 1, cntrsPlus,cntrs,ONES,sy, true);
evalRbfDer(3, 1, cntrsPlus,cntrs,ONES,sz, true);
- for (int i=0;i<numNodes ; ++i){
+ for (int i = 0; i < numNodes; ++i){
normFactor = sqrt(sx(i,0)*sx(i,0)+sy(i,0)*sy(i,0)+sz(i,0)*sz(i,0));
sx(i,0) = sx(i,0)/normFactor;
sy(i,0) = sy(i,0)/normFactor;
sz(i,0) = sz(i,0)/normFactor;
norms(i,0) = sx(i,0);norms(i,1) = sy(i,0);norms(i,2) = sz(i,0);
}
-
- for (int i=0;i<numNodes ; ++i){
- for (int j=0;j<3 ; ++j){
+
+ for (int i = 0; i < numNodes; ++i){
+ for (int j = 0; j < 3; ++j){
level_set_nodes(i,j) = cntrs(i,j);
level_set_nodes(i+numNodes,j) = cntrs(i,j)-delta*norms(i,j);
- level_set_nodes(i+2*numNodes,j) = cntrs(i,j)+delta*norms(i,j);
+ level_set_nodes(i+2*numNodes,j) = cntrs(i,j)+delta*norms(i,j);
}
level_set_funvals(i,0) = 0.0;
level_set_funvals(i+numNodes,0) = -1;
@@ -445,29 +469,31 @@ void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs,
}
}
-
-gLevelsetPoints::gLevelsetPoints(fullMatrix<double> ¢ers, int tag) : gLevelsetPrimitive(tag) {
+gLevelsetPoints::gLevelsetPoints(fullMatrix<double> ¢ers, int tag)
+ : gLevelsetPrimitive(tag)
+{
int nbNodes = 3*centers.size1();
setup_level_set(centers, points, surf);
printNodes(points, surf);
-
+
//build invA matrix for 3*n points
int indexRBF = 1;
matAInv.resize(nbNodes, nbNodes);
matAInv = generateRbfMat(0, indexRBF, points,points);
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;
}
-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");
SPoint3 sp(x,y,z);
@@ -483,35 +509,36 @@ double gLevelsetPoints::operator()(const double x, const double y, const double
// xyz_eval(0,2) = z;
// evalRbfDer(0, 1, points, xyz_eval, surf, surf_eval);
// 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);
- 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,1) = vert[i]->y();
xyz_eval(i,2) = vert[i]->z();
}
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);
}
}
/*
- assume a quadric
+ assume a quadric
x^T A x + b^T x + c = 0
- typically, we add a rotation x -> R x
+ typically, we add a rotation x -> R x
x^T (R^T A R) x + b^T R x + c = 0
and a translation x -> x - t
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++){
B[i] = lv.B[i];
for(int j = 0; j < 3; j++)
@@ -519,7 +546,9 @@ gLevelsetQuadric::gLevelsetQuadric(const gLevelsetQuadric &lv) : gLevelsetPrimit
}
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++){
res[i] = 0.;
for(int j = 0; j < 3; j++){
@@ -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){
- 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::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.);
}
-void gLevelsetQuadric::translate(const double transl[3]){
+void gLevelsetQuadric::translate(const double transl[3])
+{
double b;
xAx(transl, b, 1.0);
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]){
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.;
for(int i = 0; i < 3; i++){
b1 += B[i] * rotate[i][0];
@@ -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
-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 n[3] = {1., 0., 0.};
double ct = 1., st = 0.;
@@ -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;
}
-void gLevelsetQuadric::computeRotationMatrix(const double dir1[3], const double dir2[3], double t[3][3]){
- 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]));
+void gLevelsetQuadric::computeRotationMatrix(const double dir1[3],
+ const double dir2[3],
+ double t[3][3])
+{
+ 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 ct = 1., st = 0.;
if(norm != 0.) {
@@ -618,7 +658,8 @@ void gLevelsetQuadric::computeRotationMatrix(const double dir1[3], const double
t[2][2] = n[2] * n[2] * (1. - ct) + ct;
}
-void gLevelsetQuadric::init(){
+void gLevelsetQuadric::init()
+{
for(int i = 0; i < 3; i++){
B[i] = 0.;
for(int j = 0; j < 3; j++) A[i][j] = 0.;
@@ -626,12 +667,16 @@ void gLevelsetQuadric::init(){
C = 0.;
}
-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
- + 2. * A[1][2] * y * z + A[2][2] * z * z + B[0] * x + B[1] * y + B[2] * z + C);
+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
+ + 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
double angle = 0.;
double r;
@@ -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
double dx,dy,xi,yi,d;
dx = x-iso_x[0];
@@ -687,16 +733,20 @@ double gLevelsetShamrock::operator() (const double x, const double y, const doub
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;
sigma = _sigma;
r0 = _r0;
- xc = _xc;
+ xc = _xc;
yc = _yc;
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 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);
@@ -718,7 +768,9 @@ double gLevelsetPopcorn::operator() (const double x, const double y, const doubl
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> variables(3);
variables[0] = "x";
@@ -727,7 +779,8 @@ gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag) : gLevelsetPrimitiv
_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);
values[0] = x;
values[1] = y;
@@ -736,31 +789,31 @@ double gLevelsetMathEval::operator() (const double x, const double y, const doub
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->load(box+std::string(".msh"));
modelGeom = new GModel();
modelGeom->load(geom);
-
+
//EMI FIXME THIS
int levels = 3;
- int refineCurvedSurfaces = 0;
// double rmax = 0.1;
// double sampling = std::min(rmax, std::min(lx, std::min(ly, lz)));
-
+
//FILLING POINTS FROM GEOMBOX
std::vector<SPoint3> points;
std::vector<SPoint3> refinePoints;
for(GModel::viter vit = modelBox->firstVertex(); vit != modelBox->lastVertex(); vit++){
- for(unsigned int k = 0; k < (*vit)->getNumMeshVertices(); k++){
+ for(unsigned int k = 0; k < (*vit)->getNumMeshVertices(); k++){
MVertex *v = (*vit)->getMeshVertex(k);
SPoint3 p(v->x(), v->y(), v->z());
points.push_back(p);
}
}
for (GModel::eiter eit = modelBox->firstEdge(); eit != modelBox->lastEdge(); eit++){
- for(unsigned int k = 0; k < (*eit)->getNumMeshVertices(); k++){
+ for(unsigned int k = 0; k < (*eit)->getNumMeshVertices(); k++){
MVertex *ve = (*eit)->getMeshVertex(k);
SPoint3 pe(ve->x(), ve->y(), ve->z());
points.push_back(pe);
@@ -769,12 +822,12 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
//FILLING POINTS FROM STL
for (GModel::fiter fit = modelGeom->firstFace(); fit != modelGeom->lastFace(); fit++){
- for(unsigned int k = 0; k < (*fit)->getNumMeshVertices(); k++){
+ for(unsigned int k = 0; k < (*fit)->getNumMeshVertices(); k++){
MVertex *vf = (*fit)->getMeshVertex(k);
SPoint3 pf(vf->x(), vf->y(), vf->z());
points.push_back(pf);
}
- for(unsigned int k = 0; k < (*fit)->getNumMeshElements(); k++){
+ for(unsigned int k = 0; k < (*fit)->getNumMeshElements(); k++){
MElement *e = (*fit)->getMeshElement(k);
if (e->getType() == TYPE_TRI){
MVertex *v1 = e->getVertex(0);
@@ -794,15 +847,15 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
if (points.size() == 0) {Msg::Fatal("No points on surfaces \n"); };
SBoundingBox3d bb;
- for(unsigned int i = 0; i < points.size(); i++) bb += points[i];
- for(unsigned int i = 0; i < refinePoints.size(); i++) bb += refinePoints[i];
+ for(unsigned int i = 0; i < points.size(); i++) bb += points[i];
+ for(unsigned int i = 0; i < refinePoints.size(); i++) bb += refinePoints[i];
//bb.scale(1.01, 1.01, 1.01);
- SVector3 range = bb.max() - bb.min();
+ SVector3 range = bb.max() - bb.min();
double minLength = std::min( range.x(), std::min(range.y(), range.z()));
double hmin = minLength / 5;
- int NX = range.x() / hmin;
- int NY = range.y() / hmin;
- int NZ = range.z() / hmin;
+ int NX = range.x() / hmin;
+ int NY = range.y() / hmin;
+ int NZ = range.z() / hmin;
if(NX < 2) NX = 2;
if(NY < 2) NY = 2;
if(NZ < 2) NZ = 2;
@@ -813,7 +866,7 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
bb.max().x(), bb.max().y(), bb.max().z());
Msg::Info(" Nx=%d Ny=%d Nz=%d", NX, NY, NZ);
- _box = new cartesianBox<double>(bb.min().x(), bb.min().y(), bb.min().z(),
+ _box = new cartesianBox<double>(bb.min().x(), bb.min().y(), bb.min().z(),
SVector3(range.x(), 0, 0),
SVector3(0, range.y(), 0),
SVector3(0, 0, range.z()),
@@ -827,7 +880,7 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
while((child = parent->getChildBox())){
//Msg::Info(" level %d ", child->getLevel());
for(unsigned int i = 0; i < refinePoints.size(); i++)
- insertActiveCells(refinePoints[i].x(), refinePoints[i].y(), refinePoints[i].z(),
+ insertActiveCells(refinePoints[i].x(), refinePoints[i].y(), refinePoints[i].z(),
rtube / pow(2., (levels - child->getLevel())), *child);
parent = child;
}
@@ -839,24 +892,25 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag)
//Msg::Info("Initializing nodal values in the cartesian grid");
_box->createNodalValues();
- //Msg::Info("Computing levelset on the cartesian grid");
+ //Msg::Info("Computing levelset on the cartesian grid");
computeLevelset(modelGeom, *_box);
//Msg::Info("Renumbering mesh vertices across levels");
_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);
return dist;
}
#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()){
PView *view = PView::list[_viewIndex];
_octree = new OctreePost(view);
@@ -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.;
double val = 1.;
_octree->searchScalar(x, y, z, &val, 0);
@@ -875,8 +930,10 @@ double gLevelsetPostView::operator () (const double x, const double y, const dou
}
#endif
-gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir, const double &R,
- int tag) : gLevelsetQuadric(tag) {
+gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir,
+ const double &R, int tag)
+ : gLevelsetQuadric(tag)
+{
A[0][0] = 1.;
A[1][1] = 1.;
C = - R * R;
@@ -885,10 +942,14 @@ gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir,
rotate(rot);
translate(pt);
}
-gLevelsetGenCylinder::gLevelsetGenCylinder (const gLevelsetGenCylinder& lv) : gLevelsetQuadric(lv){}
-gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, const double &a,
- const double &b, const double &c, int tag) : gLevelsetQuadric(tag) {
+gLevelsetGenCylinder::gLevelsetGenCylinder (const gLevelsetGenCylinder& lv)
+ : gLevelsetQuadric(lv){}
+
+gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, const double &a,
+ const double &b, const double &c, int tag)
+ : gLevelsetQuadric(tag)
+{
A[0][0] = 1. / (a * a);
A[1][1] = 1. / (b * b);
A[2][2] = 1. / (c * c);
@@ -898,9 +959,13 @@ gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, cons
rotate(rot);
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[1][1] = 1.;
A[2][2] = -tan(angle) * tan(angle);
@@ -910,10 +975,15 @@ gLevelsetCone::gLevelsetCone(const double *pt, const double *dir, const double &
translate(pt);
}
-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) {
+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)
+{
A[0][0] = x2;
A[1][1] = y2;
A[2][2] = z2;
@@ -925,30 +995,37 @@ gLevelsetGeneralQuadric::gLevelsetGeneralQuadric(const double *pt, const double
translate(pt);
}
-gLevelsetGeneralQuadric::gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& lv) : gLevelsetQuadric(lv)
-{}
+gLevelsetGeneralQuadric::gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& lv)
+ : gLevelsetQuadric(lv){}
gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv) : gLevelset(lv)
{
std::vector<gLevelset *> _children=lv.getChildren();
unsigned siz = _children.size();
children.resize(siz);
- for(unsigned i = 0; i < siz; ++i)
+ for(unsigned i = 0; i < siz; ++i)
children[i] = _children[i]->clone();
}
-gLevelsetImproved::gLevelsetImproved(const gLevelsetImproved &lv) : gLevelset(lv){
+
+gLevelsetImproved::gLevelsetImproved(const gLevelsetImproved &lv)
+ : gLevelset(lv)
+{
Ls = lv.Ls->clone();
}
-gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *dir2, const double *dir3,
- const double &a, const double &b, const double &c, int tag) : gLevelsetImproved() {
+gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *dir2,
+ 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 dir2m[3] = {-dir2[0], -dir2[1], -dir2[2]};
double dir3m[3] = {-dir3[0], -dir3[1], -dir3[2]};
double n1[3]; norm(dir1, n1);
double n2[3]; norm(dir2, n2);
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]};
std::vector<gLevelset *> p;
p.push_back(new gLevelsetPlane(pt2, dir3, tag++));
@@ -960,13 +1037,18 @@ gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *d
Ls = new gLevelsetIntersection(p);
}
-gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *pt3, const double *pt4,
- const double *pt5, const double *pt6, const double *pt7, const double *pt8, int tag) : 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))
+gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *pt3,
+ const double *pt4, const double *pt5, const double *pt6,
+ const double *pt7, const double *pt8, int tag)
+ : 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",
- 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));
+ 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));
std::vector<gLevelset *> p;
p.push_back(new gLevelsetPlane(pt5, pt6, pt8, 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
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 dir1[3] = {dir[0], dir[1], dir[2]};
double dir2[3] = {-dir1[0], -dir1[1], -dir1[2]};
@@ -992,7 +1079,10 @@ gLevelsetCylinder::gLevelsetCylinder(const std::vector<double> &pt, const std::v
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 n[3]; norm(dir, n);
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
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 n[3]; norm(dir, n);
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
p2.push_back(new gLevelsetGenCylinder(pt, dir, r, tag));
Ls = new gLevelsetCut(p2);
}
-gLevelsetCylinder::gLevelsetCylinder(const gLevelsetCylinder &lv) : gLevelsetImproved(lv){}
-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() {
+gLevelsetCylinder::gLevelsetCylinder(const gLevelsetCylinder &lv)
+ : gLevelsetImproved(lv){}
+
+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 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 pt2[3] = {pt[0] + n1[0] * E - n2[0] * H2 / 2., pt[1] + n1[1] * E - n2[1] * H2 / 2.,
+ double pt1[3] = {pt[0] - n2[0] * H1 / 2., pt[1] - n2[1] * H1 / 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.};
double dir3[3]; cross(pt1, pt2, pt, dir3);
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.};
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};
@@ -1051,4 +1155,5 @@ gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1, const dou
Ls = new gLevelsetCut(p2);
}
-gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv) : gLevelsetImproved(lv){}
+gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv)
+ : gLevelsetImproved(lv){}