-/*!
- * This method works on a 3D curve linear mesh that is to say (meshDim==1 and spaceDim==3).
- * If it is not the case an exception will be thrown.
- * This method is non const because the coordinate of \a this can be appended with some new points issued from
- * intersection of plane defined by ('origin','vec').
- * This method has one in/out parameter : 'cut3DCurve'.
- * Param 'cut3DCurve' is expected to be of size 'this->getNumberOfCells()'. For each i in [0,'this->getNumberOfCells()')
- * if cut3DCurve[i]==-2, it means that for cell #i in \a this nothing has been detected previously.
- * if cut3DCurve[i]==-1, it means that cell#i has been already detected to be fully part of plane defined by ('origin','vec').
- * This method will throw an exception if \a this contains a non linear segment.
- */
-void MEDCouplingUMesh::split3DCurveWithPlane(const double *origin, const double *vec, double eps, std::vector<int>& cut3DCurve)
-{
- checkFullyDefined();
- if(getMeshDimension()!=1 || getSpaceDimension()!=3)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane works on umeshes with meshdim equal to 1 and spaceDim equal to 3 !");
- int ncells=getNumberOfCells();
- int nnodes=getNumberOfNodes();
- double vec2[3],vec3[3],vec4[3];
- double normm=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
- if(normm<1e-6)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : parameter 'vec' should have a norm2 greater than 1e-6 !");
- vec2[0]=vec[0]/normm; vec2[1]=vec[1]/normm; vec2[2]=vec[2]/normm;
- const int *conn=_nodal_connec->getConstPointer();
- const int *connI=_nodal_connec_index->getConstPointer();
- const double *coo=_coords->getConstPointer();
- std::vector<double> addCoo;
- for(int i=0;i<ncells;i++)
- {
- if(conn[connI[i]]==(int)INTERP_KERNEL::NORM_SEG2)
- {
- if(cut3DCurve[i]==-2)
- {
- int st=conn[connI[i]+1],endd=conn[connI[i]+2];
- vec3[0]=coo[3*endd]-coo[3*st]; vec3[1]=coo[3*endd+1]-coo[3*st+1]; vec3[2]=coo[3*endd+2]-coo[3*st+2];
- double normm2=sqrt(vec3[0]*vec3[0]+vec3[1]*vec3[1]+vec3[2]*vec3[2]);
- double colin=std::abs((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2])/normm2);
- if(colin>eps)//if colin<=eps -> current SEG2 is colinear to the input plane
- {
- const double *st2=coo+3*st;
- vec4[0]=st2[0]-origin[0]; vec4[1]=st2[1]-origin[1]; vec4[2]=st2[2]-origin[2];
- double pos=-(vec4[0]*vec2[0]+vec4[1]*vec2[1]+vec4[2]*vec2[2])/((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2]));
- if(pos>eps && pos<1-eps)
- {
- int nNode=((int)addCoo.size())/3;
- vec4[0]=st2[0]+pos*vec3[0]; vec4[1]=st2[1]+pos*vec3[1]; vec4[2]=st2[2]+pos*vec3[2];
- addCoo.insert(addCoo.end(),vec4,vec4+3);
- cut3DCurve[i]=nnodes+nNode;
- }
- }
- }
- }
- else
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : this method is only available for linear cell (NORM_SEG2) !");
- }
- if(!addCoo.empty())
- {
- int newNbOfNodes=nnodes+((int)addCoo.size())/3;
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo2=DataArrayDouble::New();
- coo2->alloc(newNbOfNodes,3);
- double *tmp=coo2->getPointer();
- tmp=std::copy(_coords->begin(),_coords->end(),tmp);
- std::copy(addCoo.begin(),addCoo.end(),tmp);
- DataArrayDouble::SetArrayIn(coo2,_coords);
- }
-}
-
-/*!
- * This method incarnates the policy 0 for MEDCouplingUMesh::buildExtrudedMesh method.
- * \param mesh1D is the input 1D mesh used for translation computation.
- * \return newCoords new coords filled by this method.
- */
-DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
-{
- int oldNbOfNodes=getNumberOfNodes();
- int nbOf1DCells=mesh1D->getNumberOfCells();
- int spaceDim=getSpaceDimension();
- DataArrayDouble *ret=DataArrayDouble::New();
- std::vector<bool> isQuads;
- int nbOfLevsInVec=isQuad?2*nbOf1DCells+1:nbOf1DCells+1;
- ret->alloc(oldNbOfNodes*nbOfLevsInVec,spaceDim);
- double *retPtr=ret->getPointer();
- const double *coords=getCoords()->getConstPointer();
- double *work=std::copy(coords,coords+spaceDim*oldNbOfNodes,retPtr);
- std::vector<int> v;
- std::vector<double> c;
- double vec[3];
- v.reserve(3);
- c.reserve(6);
- for(int i=0;i<nbOf1DCells;i++)
- {
- v.resize(0);
- mesh1D->getNodeIdsOfCell(i,v);
- c.resize(0);
- mesh1D->getCoordinatesOfNode(v[isQuad?2:1],c);
- mesh1D->getCoordinatesOfNode(v[0],c);
- std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
- for(int j=0;j<oldNbOfNodes;j++)
- work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
- if(isQuad)
- {
- c.resize(0);
- mesh1D->getCoordinatesOfNode(v[1],c);
- mesh1D->getCoordinatesOfNode(v[0],c);
- std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
- for(int j=0;j<oldNbOfNodes;j++)
- work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
- }
- }
- ret->copyStringInfoFrom(*getCoords());
- return ret;
-}
-
-/*!
- * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
- * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
- * \return newCoords new coords filled by this method.
- */
-DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
-{
- if(mesh1D->getSpaceDimension()==2)
- return fillExtCoordsUsingTranslAndAutoRotation2D(mesh1D,isQuad);
- if(mesh1D->getSpaceDimension()==3)
- return fillExtCoordsUsingTranslAndAutoRotation3D(mesh1D,isQuad);
- throw INTERP_KERNEL::Exception("Not implemented rotation and translation alg. for spacedim other than 2 and 3 !");
-}
-
-/*!
- * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
- * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
- * \return newCoords new coords filled by this method.
- */
-DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
-{
- if(isQuad)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : not implemented for quadratic cells !");
- int oldNbOfNodes=getNumberOfNodes();
- int nbOf1DCells=mesh1D->getNumberOfCells();
- if(nbOf1DCells<2)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
- int nbOfLevsInVec=nbOf1DCells+1;
- ret->alloc(oldNbOfNodes*nbOfLevsInVec,2);
- double *retPtr=ret->getPointer();
- retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> tmp2=getCoords()->deepCpy();
- tmp->setCoords(tmp2);
- const double *coo1D=mesh1D->getCoords()->getConstPointer();
- const int *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
- const int *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
- for(int i=1;i<nbOfLevsInVec;i++)
- {
- const double *begin=coo1D+2*conn1D[connI1D[i-1]+1];
- const double *end=coo1D+2*conn1D[connI1D[i-1]+2];
- const double *third=i+1<nbOfLevsInVec?coo1D+2*conn1D[connI1D[i]+2]:coo1D+2*conn1D[connI1D[i-2]+1];
- const double vec[2]={end[0]-begin[0],end[1]-begin[1]};
- tmp->translate(vec);
- double tmp3[2],radius,alpha,alpha0;
- const double *p0=i+1<nbOfLevsInVec?begin:third;
- const double *p1=i+1<nbOfLevsInVec?end:begin;
- const double *p2=i+1<nbOfLevsInVec?third:end;
- INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0,p1,p2,tmp3,radius,alpha,alpha0);
- double cosangle=i+1<nbOfLevsInVec?(p0[0]-tmp3[0])*(p1[0]-tmp3[0])+(p0[1]-tmp3[1])*(p1[1]-tmp3[1]):(p2[0]-tmp3[0])*(p1[0]-tmp3[0])+(p2[1]-tmp3[1])*(p1[1]-tmp3[1]);
- double angle=acos(cosangle/(radius*radius));
- tmp->rotate(end,0,angle);
- retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
- }
- return ret.retn();
-}
-
-/*!
- * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
- * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
- * \return newCoords new coords filled by this method.
- */
-DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
-{
- if(isQuad)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : not implemented for quadratic cells !");
- int oldNbOfNodes=getNumberOfNodes();
- int nbOf1DCells=mesh1D->getNumberOfCells();
- if(nbOf1DCells<2)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
- int nbOfLevsInVec=nbOf1DCells+1;
- ret->alloc(oldNbOfNodes*nbOfLevsInVec,3);
- double *retPtr=ret->getPointer();
- retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> tmp2=getCoords()->deepCpy();
- tmp->setCoords(tmp2);
- const double *coo1D=mesh1D->getCoords()->getConstPointer();
- const int *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
- const int *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
- for(int i=1;i<nbOfLevsInVec;i++)
- {
- const double *begin=coo1D+3*conn1D[connI1D[i-1]+1];
- const double *end=coo1D+3*conn1D[connI1D[i-1]+2];
- const double *third=i+1<nbOfLevsInVec?coo1D+3*conn1D[connI1D[i]+2]:coo1D+3*conn1D[connI1D[i-2]+1];
- const double vec[3]={end[0]-begin[0],end[1]-begin[1],end[2]-begin[2]};
- tmp->translate(vec);
- double tmp3[2],radius,alpha,alpha0;
- const double *p0=i+1<nbOfLevsInVec?begin:third;
- const double *p1=i+1<nbOfLevsInVec?end:begin;
- const double *p2=i+1<nbOfLevsInVec?third:end;
- double vecPlane[3]={
- (p1[1]-p0[1])*(p2[2]-p1[2])-(p1[2]-p0[2])*(p2[1]-p1[1]),
- (p1[2]-p0[2])*(p2[0]-p1[0])-(p1[0]-p0[0])*(p2[2]-p1[2]),
- (p1[0]-p0[0])*(p2[1]-p1[1])-(p1[1]-p0[1])*(p2[0]-p1[0]),
- };
- double norm=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]+vecPlane[2]*vecPlane[2]);
- if(norm>1.e-7)
- {
- vecPlane[0]/=norm; vecPlane[1]/=norm; vecPlane[2]/=norm;
- double norm2=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]);
- double vec2[2]={vecPlane[1]/norm2,-vecPlane[0]/norm2};
- double s2=norm2;
- double c2=cos(asin(s2));
- double m[3][3]={
- {vec2[0]*vec2[0]*(1-c2)+c2, vec2[0]*vec2[1]*(1-c2), vec2[1]*s2},
- {vec2[0]*vec2[1]*(1-c2), vec2[1]*vec2[1]*(1-c2)+c2, -vec2[0]*s2},
- {-vec2[1]*s2, vec2[0]*s2, c2}
- };
- double p0r[3]={m[0][0]*p0[0]+m[0][1]*p0[1]+m[0][2]*p0[2], m[1][0]*p0[0]+m[1][1]*p0[1]+m[1][2]*p0[2], m[2][0]*p0[0]+m[2][1]*p0[1]+m[2][2]*p0[2]};
- double p1r[3]={m[0][0]*p1[0]+m[0][1]*p1[1]+m[0][2]*p1[2], m[1][0]*p1[0]+m[1][1]*p1[1]+m[1][2]*p1[2], m[2][0]*p1[0]+m[2][1]*p1[1]+m[2][2]*p1[2]};
- double p2r[3]={m[0][0]*p2[0]+m[0][1]*p2[1]+m[0][2]*p2[2], m[1][0]*p2[0]+m[1][1]*p2[1]+m[1][2]*p2[2], m[2][0]*p2[0]+m[2][1]*p2[1]+m[2][2]*p2[2]};
- INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0r,p1r,p2r,tmp3,radius,alpha,alpha0);
- double cosangle=i+1<nbOfLevsInVec?(p0r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p0r[1]-tmp3[1])*(p1r[1]-tmp3[1]):(p2r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p2r[1]-tmp3[1])*(p1r[1]-tmp3[1]);
- double angle=acos(cosangle/(radius*radius));
- tmp->rotate(end,vecPlane,angle);
-
- }
- retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
- }
- return ret.retn();
-}
-
-/*!
- * This method is private because not easy to use for end user. This method is const contrary to
- * MEDCouplingUMesh::buildExtrudedMesh method because this->_coords are expected to contain
- * the coords sorted slice by slice.
- * \param isQuad specifies presence of quadratic cells.
- */
-MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMeshFromThisLowLev(int nbOfNodesOf1Lev, bool isQuad) const
-{
- int nbOf1DCells=getNumberOfNodes()/nbOfNodesOf1Lev-1;
- int nbOf2DCells=getNumberOfCells();
- int nbOf3DCells=nbOf2DCells*nbOf1DCells;
- MEDCouplingUMesh *ret=MEDCouplingUMesh::New("Extruded",getMeshDimension()+1);
- const int *conn=_nodal_connec->getConstPointer();
- const int *connI=_nodal_connec_index->getConstPointer();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn=DataArrayInt::New();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConnI=DataArrayInt::New();
- newConnI->alloc(nbOf3DCells+1,1);
- int *newConnIPtr=newConnI->getPointer();
- *newConnIPtr++=0;
- std::vector<int> newc;
- for(int j=0;j<nbOf2DCells;j++)
- {
- AppendExtrudedCell(conn+connI[j],conn+connI[j+1],nbOfNodesOf1Lev,isQuad,newc);
- *newConnIPtr++=(int)newc.size();
- }
- newConn->alloc((int)(newc.size())*nbOf1DCells,1);
- int *newConnPtr=newConn->getPointer();
- int deltaPerLev=isQuad?2*nbOfNodesOf1Lev:nbOfNodesOf1Lev;
- newConnIPtr=newConnI->getPointer();
- for(int iz=0;iz<nbOf1DCells;iz++)
- {
- if(iz!=0)
- std::transform(newConnIPtr+1,newConnIPtr+1+nbOf2DCells,newConnIPtr+1+iz*nbOf2DCells,std::bind2nd(std::plus<int>(),newConnIPtr[iz*nbOf2DCells]));
- for(std::vector<int>::const_iterator iter=newc.begin();iter!=newc.end();iter++,newConnPtr++)
- {
- int icell=(int)(iter-newc.begin());
- if(std::find(newConnIPtr,newConnIPtr+nbOf2DCells,icell)==newConnIPtr+nbOf2DCells)
- {
- if(*iter!=-1)
- *newConnPtr=(*iter)+iz*deltaPerLev;
- else
- *newConnPtr=-1;
- }
- else
- *newConnPtr=(*iter);
- }
- }
- ret->setConnectivity(newConn,newConnI,true);
- ret->setCoords(getCoords());
- return ret;
-}