Salome HOME
Contrôle des calculs de maillage
[modules/smesh.git] / src / Tools / blocFissure / gmu / meshBlocPart.py
1 # -*- coding: utf-8 -*-
2 # Copyright (C) 2014-2020  EDF R&D
3 #
4 # This library is free software; you can redistribute it and/or
5 # modify it under the terms of the GNU Lesser General Public
6 # License as published by the Free Software Foundation; either
7 # version 2.1 of the License, or (at your option) any later version.
8 #
9 # This library is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12 # Lesser General Public License for more details.
13 #
14 # You should have received a copy of the GNU Lesser General Public
15 # License along with this library; if not, write to the Free Software
16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
17 #
18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 #
20
21 import logging
22 from .geomsmesh import geompy
23 from .geomsmesh import smesh
24 from salome.smesh import smeshBuilder
25 import SMESH
26 from salome.StdMeshers import StdMeshersBuilder
27 from .putName import putName
28
29 # -----------------------------------------------------------------------------
30 # --- maillage du bloc partitionne
31
32 def meshBlocPart(blocPartition, faceFissure, tore, centres, edges, diams, circles, faces, \
33                 gencnt, facefissoutore, edgeext, facesExternes, facesExtBloc, facesExtElli, \
34                 aretesInternes, internalBoundary, ellipsoidep, sharedFaces, sharedEdges, edgesBords, \
35                 nbsegExt, nbsegGen, nbsegRad, scaleRad, reverses, reverext, nbsegCercle, nbsegFis, dmoyen, lensegEllipsoide):
36   """Maillage du bloc partitionné
37
38   TODO: a completer
39   """
40   logging.info('start')
41
42   # --- edges de bord à respecter
43
44   aFilterManager = smesh.CreateFilterManager()
45   nbAdded, internalBoundary, _NoneGroup = internalBoundary.MakeBoundaryElements( SMESH.BND_1DFROM2D, '', '', 0, [  ])
46   criteres = list()
47   unCritere = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_FreeBorders,SMESH.FT_Undefined,0)
48   criteres.append(unCritere)
49   filtre = smesh.GetFilterFromCriteria(criteres)
50   bordsLibres = internalBoundary.MakeGroupByFilter( 'bords', filtre )
51   smesh.SetName(bordsLibres, 'bordsLibres')
52
53   # --- maillage bloc
54
55   bloc1 = smesh.Mesh(blocPartition)
56
57   for i, sharedFaces_i in enumerate(sharedFaces):
58     algo2d = bloc1.Triangle(algo=smeshBuilder.NETGEN, geom=sharedFaces_i)
59     hypo2d = algo2d.Parameters(which=smesh.SIMPLE)
60     hypo2d.SetLocalLength(lensegEllipsoide)
61     hypo2d.LengthFromEdges()
62     hypo2d.SetAllowQuadrangles(0)
63     putName(algo2d.GetSubMesh(), "sharedFaces", i)
64     putName(algo2d, "algo2d_sharedFaces", i)
65     putName(hypo2d, "hypo2d_sharedFaces", i)
66
67   for i, sharedEdges_i in enumerate(sharedEdges):
68     algo1d = bloc1.Segment(geom=sharedEdges_i)
69     hypo1d = algo1d.LocalLength(lensegEllipsoide)
70     putName(algo1d.GetSubMesh(), "sharedEdges", i)
71     putName(algo1d, "algo1d_sharedEdges", i)
72     putName(hypo1d, "hypo1d_sharedEdges", i)
73
74   declareAlgoEllipsoideFirst = False
75   if declareAlgoEllipsoideFirst:
76     algo3d = bloc1.Tetrahedron(algo=smeshBuilder.NETGEN,geom=ellipsoidep)
77     hypo3d = algo3d.MaxElementVolume(1000.0)
78     putName(algo3d.GetSubMesh(), "ellipsoide")
79     putName(algo3d, "algo3d_ellipsoide")
80     putName(hypo3d, "hypo3d_ellipsoide")
81
82   algo3d = bloc1.Prism(geom=tore)
83   algo2d = bloc1.Quadrangle(geom=tore)
84   algo1d = bloc1.Segment(geom=tore)
85   hypo1d = algo1d.NumberOfSegments(nbsegGen)
86   putName(algo3d.GetSubMesh(), "tore")
87   putName(algo3d, "algo3d_tore")
88   putName(algo2d, "algo2d_tore")
89   putName(algo1d, "algo1d_tore")
90   putName(hypo1d, "hypo1d_tore")
91
92   for i, faces_i in enumerate(faces):
93     algo2d = bloc1.Quadrangle(geom=faces_i)
94     hypo2d = smesh.CreateHypothesis('QuadrangleParams')
95     hypo2d.SetTriaVertex( geompy.GetSubShapeID(blocPartition,centres[i]) )
96     hypo2d.SetQuadType( StdMeshersBuilder.QUAD_STANDARD )
97     status = bloc1.AddHypothesis(hypo2d,faces_i)
98     putName(algo2d.GetSubMesh(), "faces", i)
99     putName(algo2d, "algo2d_faces", i)
100     putName(hypo2d, "hypo2d_faces", i)
101
102   for i, edges_i in enumerate(edges):
103     algo1d = bloc1.Segment(geom=edges_i)
104     if reverses[i] > 0:
105       hypo1d = algo1d.NumberOfSegments(nbsegRad, scaleRad,[ geompy.GetSubShapeID(blocPartition,edges_i) ])
106     else:
107       hypo1d = algo1d.NumberOfSegments(nbsegRad, scaleRad,[  ])
108     putName(algo1d.GetSubMesh(), "edges", i)
109     putName(algo1d, "algo1d_edges", i)
110     putName(hypo1d, "hypo1d_edges", i)
111
112   for i, circles_i in enumerate(circles):
113     algo1d = bloc1.Segment(geom=circles_i)
114     hypo1d = algo1d.NumberOfSegments(nbsegCercle)
115     putName(algo1d.GetSubMesh(), "circles", i)
116     putName(algo1d, "algo1d_circles", i)
117     putName(hypo1d, "hypo1d_circles", i)
118
119   if len(edgeext) == 1:
120     densite = int(round(nbsegFis/2))
121     algo1d = bloc1.Segment(geom=edgeext[0])
122     hypo1d = algo1d.NumberOfSegments(nbsegFis)
123     hypo1d.SetDistrType( 2 )
124     hypo1d.SetConversionMode( 1 )
125     hypo1d.SetTableFunction( [ 0, densite, 0.4, 1, 0.6, 1, 1, densite ] )
126     putName(algo1d.GetSubMesh(), "edgeext")
127     putName(algo1d, "algo1d_edgeext")
128     putName(hypo1d, "hypo1d_edgeext")
129   else:
130     longTotal = 0
131     longEdgeExts = list()
132     for i, edgeext_i in enumerate(edgeext):
133       props = geompy.BasicProperties(edgeext_i)
134       longEdgeExts.append(props[0])
135       longTotal += props[0]
136     for i, edgeext_i in enumerate(edgeext):
137       local = longTotal/nbsegFis
138       nbLocal = int(round(nbsegFis*longEdgeExts[i]/longTotal))
139       densite = int(round(nbLocal/2))
140       algo1d = bloc1.Segment(geom=edgeext_i)
141       hypo1d = algo1d.NumberOfSegments(nbLocal)
142       hypo1d.SetDistrType( 2 )
143       hypo1d.SetConversionMode( 1 )
144       hypo1d.SetTableFunction( [ 0, densite, 0.8, 1, 1, 1 ] )
145       if reverext[i]:
146         hypo1d.SetReversedEdges([ geompy.GetSubShapeID(blocPartition, edgeext_i) ])
147       putName(algo1d.GetSubMesh(), "edgeext", i)
148       putName(algo1d, "algo1d_edgeext", i)
149       putName(hypo1d, "hypo1d_edgeext", i)
150
151   algo2d = bloc1.Triangle(algo=smeshBuilder.NETGEN_2D, geom=facefissoutore)
152   hypo2d = algo2d.LengthFromEdges()
153   putName(algo2d.GetSubMesh(), "facefissoutore")
154   putName(algo2d, "algo2d_facefissoutore")
155   putName(hypo2d, "hypo2d_facefissoutore")
156
157
158   maxElemArea = 0.5*dmoyen*dmoyen
159   logging.debug("dmoyen %s, maxElemArea %s", dmoyen, maxElemArea)
160
161   for i, facesExternes_i in enumerate(facesExternes):
162     algo2d = bloc1.Triangle(algo=smeshBuilder.NETGEN_2D, geom=facesExternes_i)
163     hypo2d = algo2d.MaxElementArea(maxElemArea)
164     if edgesBords is None:
165       algo1d = bloc1.Segment(geom=facesExternes_i)
166       hypo1d = algo1d.NumberOfSegments(1)
167     putName(algo2d.GetSubMesh(), "facesExternes", i)
168     putName(algo2d, "algo2d_facesExternes", i)
169     putName(hypo2d, "hypo2d_facesExternes", i)
170     if edgesBords is None:
171       putName(algo1d, "algo1d_facesExternes", i)
172       putName(hypo1d, "hypo1d_facesExternes", i)
173
174   for i, aretesInternes_i in enumerate(aretesInternes):
175     algo1d = bloc1.Segment(geom=aretesInternes_i)
176     hypo1d = algo1d.NumberOfSegments(nbsegExt)
177     putName(algo1d.GetSubMesh(), "aretesInternes", i)
178     putName(algo1d, "algo1d_aretesInternes", i)
179     putName(hypo1d, "hypo1d_aretesInternes", i)
180
181   if edgesBords is not None:
182     algo1d = bloc1.UseExisting1DElements(geom=edgesBords)
183     hypo1d = algo1d.SourceEdges([ bordsLibres ],0,0)
184     putName(algo1d.GetSubMesh(), "bordsLibres")
185     putName(algo1d, "algo1d_bordsLibres")
186     putName(hypo1d, "hypo1d_bordsLibres")
187
188   #isDone = bloc1.Compute()
189
190   if not declareAlgoEllipsoideFirst:
191     algo3d = bloc1.Tetrahedron(algo=smeshBuilder.NETGEN,geom=ellipsoidep)
192     hypo3d = algo3d.MaxElementVolume(1000.0)
193     putName(algo3d.GetSubMesh(), "ellipsoide")
194     putName(algo3d, "algo3d_ellipsoide")
195     putName(hypo3d, "hypo3d_ellipsoide")
196
197   is_done = bloc1.Compute()
198   text = "bloc1.Compute"
199   if is_done:
200     logging.info(text+" OK")
201   else:
202     text = "Erreur au calcul du maillage.\n" + text
203     logging.info(text)
204     raise Exception(text)
205
206   nbRemoved = bloc1.RemoveOrphanNodes()
207
208   faceFissure1 = bloc1.GroupOnGeom(faceFissure,'FACE1',SMESH.FACE)
209   noeudsFondFissure = bloc1.GroupOnGeom(gencnt,'nfondfis',SMESH.NODE)
210
211   groups_faceCommuneEllipsoideBloc = list()
212   for i, sharedFaces_i in enumerate(sharedFaces):
213     name = "faceCommuneEllipsoideBloc_%d"%i
214     groups_faceCommuneEllipsoideBloc.append(bloc1.GroupOnGeom(sharedFaces_i, name, SMESH.FACE))
215   groups_faceExterneBloc = list()
216   for i, facesExtBloc_i in enumerate(facesExtBloc):
217     name = "faceExterneBloc_%d"%i
218     groups_faceExterneBloc.append(bloc1.GroupOnGeom(facesExtBloc_i, name, SMESH.FACE))
219
220   skinBlocMeshes = list()
221   for i, groups_faceCommuneEllipsoideBloc_i in enumerate(groups_faceCommuneEllipsoideBloc):
222     name = "faceCommuneEllipsoideBloc_%d"%i
223     skinBlocMeshes.append(smesh.CopyMesh(groups_faceCommuneEllipsoideBloc_i, name, 0, 0))
224   for i, groups_faceExterneBloc_i in enumerate(groups_faceExterneBloc):
225     name = "faceExterneBloc_%d"%i
226     skinBlocMeshes.append(smesh.CopyMesh(groups_faceExterneBloc_i, name, 0, 0))
227
228   meshesBloc = [internalBoundary.GetMesh()]
229   for i, skinBlocMeshes_i in enumerate(skinBlocMeshes):
230     meshesBloc.append(skinBlocMeshes_i.GetMesh())
231   blocMesh = smesh.Concatenate(meshesBloc, 1, 1, 1e-05,False)
232
233   algo3d = blocMesh.Tetrahedron(algo=smeshBuilder.NETGEN)
234   hypo3d = algo3d.MaxElementVolume(1000.0)
235   putName(algo3d.GetSubMesh(), "bloc")
236   putName(algo3d, "algo3d_bloc")
237   putName(hypo3d, "hypo3d_bloc")
238
239   is_done = blocMesh.Compute()
240   text = "blocMesh.Compute"
241   if is_done:
242     logging.info(text+" OK")
243   else:
244     text = "Erreur au calcul du maillage.\n" + text
245     logging.info(text)
246     raise Exception(text)
247
248   blocComplet = smesh.Concatenate([bloc1.GetMesh(), blocMesh.GetMesh()], 1, 1, 1e-05,False)
249
250   return bloc1, blocComplet