for(int i=0;i<nbOfCells;i++,nodalI++,retPtr+=4)
{
double matrix[16]={0,0,0,1,0,0,0,1,0,0,0,1,1,1,1,0},matrix2[16];
- if(nodalI[1]-nodalI[0]>=3)
+ if(nodalI[1]-nodalI[0]>=4)
{
+ double aa[3]={coor[nodal[nodalI[0]+1+1]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0],
+ coor[nodal[nodalI[0]+1+1]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1],
+ coor[nodal[nodalI[0]+1+1]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]}
+ ,bb[3]={coor[nodal[nodalI[0]+1+2]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0],
+ coor[nodal[nodalI[0]+1+2]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1],
+ coor[nodal[nodalI[0]+1+2]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]};
+ double cc[3]={aa[1]*bb[2]-aa[2]*bb[1],aa[2]*bb[0]-aa[0]*bb[2],aa[0]*bb[1]-aa[1]*bb[0]};
for(int j=0;j<3;j++)
{
int nodeId(nodal[nodalI[0]+1+j]);
throw INTERP_KERNEL::Exception(oss.str());
}
}
+ if(sqrt(cc[0]*cc[0]+cc[1]*cc[1]+cc[2]*cc[2])>1e-7)
+ {
+ INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
+ retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
+ }
+ else
+ {
+ if(nodalI[1]-nodalI[0]==4)
+ {
+ std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : cell" << i << " : Presence of The 3 colinear points !";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+ //
+ double dd[3]={0.,0.,0.};
+ for(int offset=nodalI[0]+1;offset<nodalI[1];offset++)
+ std::transform(coor+3*nodal[offset],coor+3*(nodal[offset]+1),dd,dd,std::plus<double>());
+ int nbOfNodesInCell(nodalI[1]-nodalI[0]-1);
+ std::transform(dd,dd+3,dd,std::bind2nd(std::multiplies<double>(),1./(double)nbOfNodesInCell));
+ std::copy(dd,dd+3,matrix+4*2);
+ INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
+ retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
+ }
}
else
{
std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : invalid 2D cell #" << i << " ! Must be constitued by more than 3 nodes !";
throw INTERP_KERNEL::Exception(oss.str());
}
- INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
- retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
}
return ret.retn();
}
return MergeVorCells2D(p,eps,true);
}
+/*!
+ * suppress additional sub points on edges
+ */
+MCAuto<MEDCouplingUMesh> SimplifyPolygon(const MEDCouplingUMesh *m, double eps)
+{
+ if(m->getNumberOfCells()!=1)
+ throw INTERP_KERNEL::Exception("SimplifyPolygon : internal error !");
+ const int *conn(m->getNodalConnectivity()->begin()),*conni(m->getNodalConnectivityIndex()->begin());
+ int nbPtsInPolygon(conni[1]-conni[0]-1);
+ const double *coo(m->getCoords()->begin());
+ std::vector<int> resConn;
+ for(int i=0;i<nbPtsInPolygon;i++)
+ {
+ int prev(conn[(i+nbPtsInPolygon-1)%nbPtsInPolygon+1]),current(conn[i%nbPtsInPolygon+1]),zeNext(conn[(i+1)%nbPtsInPolygon+1]);
+ double a[3]={
+ coo[3*prev+0]-coo[3*current+0],
+ coo[3*prev+1]-coo[3*current+1],
+ coo[3*prev+2]-coo[3*current+2],
+ },b[3]={
+ coo[3*current+0]-coo[3*zeNext+0],
+ coo[3*current+1]-coo[3*zeNext+1],
+ coo[3*current+2]-coo[3*zeNext+2],
+ };
+ double c[3]={a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]};
+ if(sqrt(c[0]*c[0]+c[1]*c[1]+c[2]*c[2])>eps)
+ resConn.push_back(current);
+ }
+ MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
+ ret->setCoords(m->getCoords());
+ ret->allocateCells();
+ ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,resConn.size(),&resConn[0]);
+ return ret;
+}
+
MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
{
std::size_t sz(vcs.size());
tmp2=tmp;
else
tmp2=MergeVorCells2D(tmp,eps,false);
+ tmp2=SimplifyPolygon(tmp2,eps);
const int *cPtr(tmp2->getNodalConnectivity()->begin()),*ciPtr(tmp2->getNodalConnectivityIndex()->begin());
conn.insert(conn.end(),cPtr+1,cPtr+ciPtr[1]);
}
vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
if(polygsToIterOn.size()<1)
throw INTERP_KERNEL::Exception("Voronoize3D : presence of a point outside the given cell !");
- std::set<int> elemsToDo(polygsToIterOn.begin(),polygsToIterOn.end()),elemsDone;
- std::size_t ii(0);
std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
- MCAuto<DataArrayInt> d(DataArrayInt::New()),dI(DataArrayInt::New()),rd(DataArrayInt::New()),rdI(DataArrayInt::New());
- MCAuto<MEDCouplingUMesh> faces(vorTess->buildDescendingConnectivity(d,dI,rd,rdI));
- //
- while(!elemsToDo.empty())
+ for(int poly=0;poly<vorTess->getNumberOfCells();poly++)
{
- int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
const double *seed(pts+3*poly);
MCAuto<MEDCouplingUMesh> tile(l0[poly]);
tile->zipCoords();
newVorCell->zipCoords();
MCAuto<MEDCouplingUMesh> modifiedCell(cells->buildPartOfMySelfSlice(0,1,1,true));
modifiedCell->zipCoords();
- l0[poly]=modifiedCell;
- if(std::find(polygsToIterOn.begin(),polygsToIterOn.end(),poly)!=polygsToIterOn.end())// we iterate on a polyhedron containg the point to add pt -> add cells sharing faces with just computed newVorCell
- {
- MCAuto<MEDCouplingUMesh> faces2;
- {
- MCAuto<DataArrayInt> d2(DataArrayInt::New()),d2I(DataArrayInt::New()),rd2(DataArrayInt::New()),rd2I(DataArrayInt::New());
- faces2=newVorCell->buildDescendingConnectivity(d2,d2I,rd2,rd2I);
- }
- MCAuto<MEDCouplingUMesh> faces3(faces2->buildPartOfMySelfSlice(1,faces2->getNumberOfCells(),1,true));// suppress internal face
- MCAuto<MEDCouplingUMesh> facesOfCurSplitPol(faces->buildPartOfMySelf(d->begin()+dI->getIJ(poly,0),d->begin()+dI->getIJ(poly+1,0),true));
- // intersection between the out faces of newVorCell and the neighbor faces of poly polyhedron -> candidates
- MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(facesOfCurSplitPol);
- MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(faces3);
- INTERP_KERNEL::Interpolation3DSurf interpolation;
- interpolation.setMinDotBtwPlane3DSurfIntersect(eps2);
- interpolation.setMaxDistance3DSurfIntersect(eps);
- interpolation.setPrecision(1e-12);
- std::vector<std::map<int,double> > matrix;
- interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix,"P0P0");
- std::set<int> zeCandidates;
- for(std::vector<std::map<int,double> >::const_iterator it2=matrix.begin();it2!=matrix.end();it2++)
- for(std::map<int,double>::const_iterator it3=(*it2).begin();it3!=(*it2).end();it3++)
- {
- int faceIdInVorTess(d->getIJ(dI->getIJ(poly,0)+(*it3).first,0));
- for(const int *it4=rd->begin()+rdI->getIJ(faceIdInVorTess,0);it4!=rd->begin()+rdI->getIJ(faceIdInVorTess+1,0);it4++)
- {
- if(*it4!=poly)
- zeCandidates.insert(*it4);
- }
- }
- std::set<int> tmp,newElementsToDo;
- std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
- std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
- elemsToDo=newElementsToDo;
- }
- //
newVorCells.push_back(newVorCell);
- ii++;
+ l0[poly]=modifiedCell;
}
l0.push_back(MergeVorCells3D(newVorCells,eps));
}