src/Plugins/DevelopedSurface/scripts/pp2.py

   1 # Copyright (C) 2017-2023  CEA/DEN, EDF R&D
   2 #
   3 # This library is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU Lesser General Public
   5 # License as published by the Free Software Foundation; either
   6 # version 2.1 of the License, or (at your option) any later version.
   7 #
   8 # This library is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 # Lesser General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU Lesser General Public
  14 # License along with this library; if not, write to the Free Software
  15 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 #
  17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 #
  19 # Author : Anthony Geay (EDF R&D)
  20
  21 from MEDLoader import *
  22 from math import pi,asin,acos,atan2
  23 import numpy as np
  24
  25 invertThetaInc=False
  26
  27 mm=MEDFileMesh.New("example2_2.med")
  28 c=mm.getCoords()
  29
  30 v=DataArrayDouble([0.70710678118617443,-0.70710678118572623,1],1,3)
  31 v/=v.magnitude()[0]
  32 vRef=DataArrayDouble([0,0,1],1,3)
  33 vRot=DataArrayDouble.CrossProduct(v,vRef)
  34 ang=asin(vRot.magnitude()[0])
  35 MEDCouplingUMesh.Rotate3DAlg([0,0,0],vRot.getValues(),ang,c)
  36
  37 c2=c[:,[0,1]]
  38 #
  39 a,b=c2.findCommonTuples(1e-5)
  40 pool=c2[a[b[:-1]]]
  41 probes=[pool[:,0].getMaxValue()[1],pool[:,1].getMaxValue()[1],pool[:,0].getMinValue()[1]]
  42 center,radius,_=pool[probes].asArcOfCircle()
  43 assert(((c2-center).magnitude()).isUniform(radius,1e-5))
  44 theta=((c2-center)/radius)
  45 theta=DataArrayDouble([atan2(y,x) for x,y in theta])
  46 zeTheta=theta-theta[0]
  47
  48 mm.write("tmp.med",2)
  49 ###########
  50 c_cyl=DataArrayDouble.fromCartToCyl(c-(center+[0.]))
  51 assert(c_cyl[:,0].isUniform(radius,1e-5))
  52 tetha0=c_cyl[0,1]
  53 c_cyl[:,1]-=tetha0
  54 c_cyl_2=c_cyl[:,1] ; c_cyl_2.abs()
  55 c_cyl[c_cyl_2.findIdsInRange(0.,1e-6),1]=0.
  56 ########
  57 m0=mm[0]
  58 c0s=m0.getCellIdsLyingOnNodes([0],False)
  59 assert(len(c0s)!=0)
  60 tmp=c_cyl[m0.getNodeIdsOfCell(c0s[0]),1].getMaxAbsValueInArray()
  61 a=tmp>0.
  62 if a^(not invertThetaInc):
  63     c_cyl[:,1]=-c_cyl[:,1]
  64     pass
  65 #
  66 c_cyl_post=c_cyl[:,1]
  67 m0.convertAllToPoly()
  68 tmp=MEDCoupling1DGTUMesh(m0) ; c=tmp.getNodalConnectivity() ; ci=tmp.getNodalConnectivityIndex()
  69 for elt in c_cyl_post:
  70     if float(elt)<0:
  71         elt[:]=float(elt)+2*pi
  72
  73 nbCells=m0.getNumberOfCells()
  74 nbNodes=m0.getNumberOfNodes()
  75 newCoo=DataArrayInt(0,1)
  76 cellsWithPb=[]
  77 newConn=[]
  78 for i in xrange(nbCells):
  79     tmp=c_cyl_post[c[ci[i]:ci[i+1]]]
  80     if tmp.getMaxValueInArray()-tmp.getMinValueInArray()>pi:
  81         cellsWithPb.append(i)
  82         # dup of low val
  83         newLocConn=[]
  84         for elt in c[ci[i]:ci[i+1]]:
  85             if float(c_cyl_post[elt])<=pi:
  86                 newLocConn.append(nbNodes+len(newCoo))
  87                 newCoo.pushBackSilent(int(elt))
  88                 pass
  89             else:
  90                 newLocConn.append(int(elt))
  91                 pass
  92             pass
  93         newConn.append(newLocConn[:])
  94         pass
  95
  96
  97 c_cyl_post=DataArrayDouble.Aggregate([c_cyl_post,c_cyl_post[newCoo]+2*pi])
  98 z_post=DataArrayDouble.Aggregate([c_cyl[:,2],c_cyl[newCoo,2]])
  99 newCoords=DataArrayDouble.Meld(radius*c_cyl_post,z_post)
 100
 101 z=DataArrayInt([len(elt) for elt in newConn])
 102 z.computeOffsetsFull()
 103 #
 104 cellsWithPbMesh=MEDCoupling1DGTUMesh("",NORM_POLYGON)
 105 cellsWithPbMesh.setNodalConnectivity(DataArrayInt(sum(newConn,[])),z)
 106 cellsWithPbMesh.setCoords(newCoords)
 107 #
 108 ko_part_ids=DataArrayInt(cellsWithPb)
 109 part_ids=ko_part_ids.buildComplement(nbCells)
 110 part=m0[part_ids]
 111 part.setCoords(newCoords)
 112 whole=MEDCouplingUMesh.MergeUMeshesOnSameCoords([part,cellsWithPbMesh.buildUnstructured()])
 113
 114 o2n=DataArrayInt.Aggregate([part_ids,ko_part_ids])
 115 whole=whole[o2n.invertArrayO2N2N2O(len(o2n))]
 116 whole.setName("mesh")
 117 WriteMesh("tmp3.med",whole,True)
 118
 119 #sk=m0.computeSkin().computeFetchedNodeIds()
 120 #sk_part=part.computeSkin().computeFetchedNodeIds()
 121 #m0[DataArrayInt(cellsWithPb)].computeFetchedNodeIds()