*/
DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception)
{
- throw INTERP_KERNEL::Exception("Not implemented yet !");
+ 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 !");
+ 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);
+ MEDCouplingUMesh *tmp=MEDCouplingUMesh::New();
+ DataArrayDouble *tmp2=getCoords()->deepCopy();
+ tmp->setCoords(tmp2);
+ tmp2->decrRef();
+ 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);
+ }
+ tmp->decrRef();
+ return ret;
}
/*!
self.assertAlmostEqual(expected2[i],f2.getIJ(0,i),12);
pass
pass
+
+ def testExtrudedMesh7(self):
+ coo1=[0.,1.,2.,3.5]
+ a=DataArrayDouble.New();
+ a.setValues(coo1,4,1);
+ b=MEDCouplingCMesh.New();
+ b.setCoordsAt(0,a);
+ c=b.buildUnstructured();
+ self.assertEqual(1,c.getSpaceDimension());
+ c.changeSpaceDimension(2);
+ #
+ d=DataArrayDouble.New();
+ d.alloc(13,1);
+ d.iota();
+ e=MEDCouplingCMesh.New();
+ e.setCoordsAt(0,d);
+ f=e.buildUnstructured();
+ g=f.getCoords().applyFunc(2,"3.5*IVec+x/6*3.14159265359*JVec");
+ h=g.fromPolarToCart();
+ f.setCoords(h);
+ i=c.buildExtrudedMeshFromThis(f,1);
+ self.assertEqual(52,i.getNumberOfNodes());
+ tmp,tmp2,tmp3=i.mergeNodes(1e-9);
+ self.assertTrue(tmp2);
+ self.assertEqual(37,tmp3);
+ i.convertDegeneratedCells();
+ vec1=[10.,0.,0.]
+ i.translate(vec1);
+ g2=h.applyFunc(3,"13.5/3.5*x*IVec+0*JVec+13.5/3.5*y*KVec");
+ f.setCoords(g2);
+ i.changeSpaceDimension(3);
+ i3=i.buildExtrudedMeshFromThis(f,1);
+ f2=i3.getMeasureField(True);
+ tmp,tmp2,tmp3=i.mergeNodes(1e-9);
+ self.assertTrue(tmp2);
+ self.assertEqual(444,tmp3);
+ expected1=[1.327751058489274, 4.2942574094314701, 13.024068164857139, 1.3069177251569044, 4.1484240761012954, 12.297505664866796, 1.270833333332571, 3.8958333333309674, 11.039062499993179, 1.2291666666659207, 3.6041666666644425, 9.585937499993932, 1.1930822748415895, 3.3515759238941376, 8.3274943351204556, 1.1722489415082769, 3.2057425905609289, 7.6009318351210622, 1.1722489415082862, 3.2057425905609884, 7.6009318351213713, 1.1930822748416161, 3.3515759238943001, 8.3274943351212727, 1.2291666666659564, 3.6041666666646734, 9.5859374999950777, 1.2708333333326081, 3.8958333333311868, 11.039062499994293, 1.3069177251569224, 4.1484240761014384, 12.297505664867627, 1.3277510584902354, 4.2942574094346071, 13.024068164866796]
+ for ii in xrange(12):
+ for jj in xrange(36):
+ self.assertAlmostEqual(expected1[jj],f2.getIJ(0,ii*36+jj),9);
+ pass
+ #
+ pass
def setUp(self):
pass
