From: ageay <ageay>
Date: Thu, 2 Dec 2010 10:42:49 +0000 (+0000)
Subject: *** empty log message ***
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=87132069be97288dc9b9a7093a19ae6a07b4caaa;p=tools%2Fmedcoupling.git

*** empty log message ***
---

diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index 9db650a54..17633fa23 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -2173,7 +2173,66 @@ DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D(con
  */
 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;
 }
 
 /*!
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index 64515da85..6633b2531 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -5046,6 +5046,49 @@ class MEDCouplingBasicsTest(unittest.TestCase):
             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