Arcane  v3.16.0.0
Documentation développeur
Chargement...
Recherche...
Aucune correspondance
UnstructuredMeshUtilities.cc
1// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2//-----------------------------------------------------------------------------
3// Copyright 2000-2025 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4// See the top-level COPYRIGHT file for details.
5// SPDX-License-Identifier: Apache-2.0
6//-----------------------------------------------------------------------------
7/*---------------------------------------------------------------------------*/
8/* UnstructuredMeshUtilities.cc (C) 2000-2025 */
9/* */
10/* Fonctions utilitaires sur un maillage. */
11/*---------------------------------------------------------------------------*/
12/*---------------------------------------------------------------------------*/
13
14#include "arcane/utils/InvalidArgumentException.h"
15#include "arcane/utils/FatalErrorException.h"
16#include "arcane/utils/ParallelFatalErrorException.h"
17#include "arcane/utils/ArgumentException.h"
18#include "arcane/utils/Array.h"
19#include "arcane/utils/Real3.h"
20#include "arcane/utils/HashTableMap.h"
21#include "arcane/utils/OStringStream.h"
22#include "arcane/utils/ScopedPtr.h"
23#include "arcane/utils/StringBuilder.h"
24#include "arcane/utils/CheckedConvert.h"
25#include "arcane/utils/PlatformUtils.h"
26#include "arcane/utils/SmallArray.h"
27#include "arcane/utils/NotImplementedException.h"
28
29#include "arcane/IMesh.h"
30#include "arcane/IMeshWriter.h"
31#include "arcane/IParallelMng.h"
32#include "arcane/MeshUtils.h"
33#include "arcane/IItemFamily.h"
34#include "arcane/ItemPrinter.h"
35#include "arcane/VariableTypes.h"
36#include "arcane/ItemPairGroup.h"
37#include "arcane/ISubDomain.h"
38#include "arcane/ServiceBuilder.h"
39#include "arcane/IParallelReplication.h"
40#include "arcane/Timer.h"
41#include "arcane/IMeshPartitioner.h"
42#include "arcane/IPrimaryMesh.h"
43#include "arcane/IMeshChecker.h"
44
45#include "arcane/mesh/UnstructuredMeshUtilities.h"
46#include "arcane/IItemFamilyNetwork.h"
47#include "arcane/ConnectivityItemVector.h"
48#include "arcane/mesh/NewItemOwnerBuilder.h"
49#include "arcane/mesh/ParticleFamily.h"
50#include "arcane/mesh/GraphDoFs.h"
51#include "arcane/mesh/BasicItemPairGroupComputeFunctor.h"
52#include "arcane/mesh/MeshNodeMerger.h"
53#include "arcane/mesh/ConnectivityNewWithDependenciesTypes.h"
54#include "arcane/mesh/ItemsOwnerBuilder.h"
55
56#include <algorithm>
57
58/*---------------------------------------------------------------------------*/
59/*---------------------------------------------------------------------------*/
60
61namespace Arcane
62{
63
64/*---------------------------------------------------------------------------*/
65/*---------------------------------------------------------------------------*/
66
67UnstructuredMeshUtilities::
68UnstructuredMeshUtilities(IMesh* mesh)
69: TraceAccessor(mesh->traceMng())
70, m_mesh(mesh)
71, m_compute_adjacency_functor(new BasicItemPairGroupComputeFunctor(traceMng()))
72{
73}
74
75/*---------------------------------------------------------------------------*/
76/*---------------------------------------------------------------------------*/
77
78UnstructuredMeshUtilities::
79~UnstructuredMeshUtilities()
80{
81 delete m_compute_adjacency_functor;
82}
83
84/*---------------------------------------------------------------------------*/
85/*---------------------------------------------------------------------------*/
86
87void UnstructuredMeshUtilities::
88changeOwnersFromCells()
89{
90 // On suppose qu'on connait les nouveaux propriétaires des mailles, qui
91 // se trouvent dans cells_owner. Il faut
92 // maintenant déterminer les nouveaux propriétaires des noeuds et
93 // des faces. En attendant d'avoir un algorithme qui équilibre mieux
94 // les messages, on applique le suivant:
95 // - chaque sous-domaine est responsable pour déterminer le nouveau
96 // propriétaire des noeuds et des faces qui lui appartiennent.
97 // - pour les noeuds, le nouveau propriétaire est le nouveau propriétaire
98 // de la maille connectée à ce noeud dont le uniqueId() est le plus petit.
99 // - pour les faces, le nouveau propriétaire est le nouveau propriétaire
100 // de la maille qui est derrière cette face s'il s'agit d'une face
101 // interne et de la maille connectée s'il s'agit d'une face frontière.
102 // - pour les noeuds duaux, le nouveau propriétaire est le nouveau propriétaire
103 // de la maille connectée à l'élément dual
104 // - pour les liaisons, le nouveau propriétaire est le nouveau propriétaire
105 // de la maille connectée au premier noeud dual, c'est-à-dire le propriétaire
106 // du premier noeud dual de la liaison
107
108 // Outil d'affectation des owners pour les items
109 mesh::NewItemOwnerBuilder owner_builder;
110
111 VariableItemInt32& nodes_owner(m_mesh->nodeFamily()->itemsNewOwner());
112 VariableItemInt32& edges_owner(m_mesh->edgeFamily()->itemsNewOwner());
113 VariableItemInt32& faces_owner(m_mesh->faceFamily()->itemsNewOwner());
114 VariableItemInt32& cells_owner(m_mesh->cellFamily()->itemsNewOwner());
115
116 // Détermine les nouveaux propriétaires des noeuds
117 {
118 ENUMERATE_NODE(i_node,m_mesh->ownNodes()){
119 const Node node = *i_node;
120 const Cell cell = owner_builder.connectedCellOfItem(node);
121#ifdef ARCANE_DEBUG_LOAD_BALANCING
122 if (nodes_owner[node]!= cells_owner[cell]){
123 info() << "New owner for node: " << ItemPrinter(node) << " cell=" << ItemPrinter(cell)
124 << " old_owner=" << nodes_owner[node]
125 << " current_cell_owner=" << cell.owner()
126 << " new_owner=" << cells_owner[cell];
127 }
128#endif /* ARCANE_DEBUG_LOAD_BALANCING */
129 nodes_owner[node] = cells_owner[cell];
130 }
131 nodes_owner.synchronize();
132 }
133
134 // Détermine les nouveaux propriétaires des arêtes
135 {
136 ENUMERATE_EDGE(i_edge,m_mesh->ownEdges()){
137 const Edge edge = *i_edge;
138 const Cell cell = owner_builder.connectedCellOfItem(edge);
139#ifdef ARCANE_DEBUG_LOAD_BALANCING
140 if (edges_owner[edge] != cells_owner[cell]) {
141 info() << "New owner for edge: " << ItemPrinter(edge) << " cell=" << ItemPrinter(cell)
142 << " old_owner=" << edges_owner[edge]
143 << " current_cell_owner=" << cell.owner()
144 << " new_owner=" << cells_owner[cell];
145 }
146#endif /* ARCANE_DEBUG_LOAD_BALANCING */
147 edges_owner[edge] = cells_owner[cell];
148 }
149 edges_owner.synchronize();
150 }
151
152 // Détermine les nouveaux propriétaires des faces
153 {
154 ENUMERATE_FACE(i_face,m_mesh->ownFaces()){
155 const Face face = *i_face;
156 const Cell cell = owner_builder.connectedCellOfItem(face);
157 faces_owner[face] = cells_owner[cell];
158 }
159 faces_owner.synchronize();
160 }
161
162
163 // Détermine les nouveaux propriétaires des particules
164 // Les particules ont le même propriétaire que celui de la maille dans
165 // laquelle elle se trouve.
166 for( IItemFamily* family : m_mesh->itemFamilies() ){
167 if (family->itemKind()!=IK_Particle)
168 continue;
169 // Positionne les nouveaux propriétaires des particle
170 VariableItemInt32& particles_owner(family->itemsNewOwner());
171 ENUMERATE_PARTICLE(i_particle,family->allItems()){
172 Particle particle = *i_particle ;
173 particles_owner[particle] = cells_owner[particle.cell()] ;
174 }
175 }
176
177 // GraphOnDoF
178 if(m_mesh->itemFamilyNetwork())
179 {
180 // Dof with Mesh Item connectivity
181 for( IItemFamily* family : m_mesh->itemFamilies() )
182 {
183 if (family->itemKind()!=IK_DoF || family->name()==mesh::GraphDoFs::linkFamilyName())
184 continue;
185 VariableItemInt32& dofs_new_owner(family->itemsNewOwner());
186 std::array<Arcane::eItemKind,5> dualitem_kinds = {IK_Cell,IK_Face,IK_Edge,IK_Node,IK_Particle} ;
187 for(auto dualitem_kind : dualitem_kinds)
188 {
189 IItemFamily* dualitem_family = dualitem_kind==IK_Particle? m_mesh->findItemFamily(dualitem_kind,mesh::ParticleFamily::defaultFamilyName(), false):
190 m_mesh->itemFamily(dualitem_kind) ;
191 if(dualitem_family)
192 {
193
194 VariableItemInt32& dualitems_new_owner(dualitem_family->itemsNewOwner());
195 auto connectivity_name = mesh::connectivityName(family,dualitem_family) ;
196 bool is_dof2dual = true ;
197 auto connectivity = m_mesh->itemFamilyNetwork()->getConnectivity(family,dualitem_family,connectivity_name) ;
198 if(!connectivity)
199 {
200 connectivity = m_mesh->itemFamilyNetwork()->getConnectivity(dualitem_family,family,connectivity_name) ;
201 is_dof2dual = false ;
202 }
203
204 if(connectivity)
205 {
206 ConnectivityItemVector accessor(connectivity);
207 if(is_dof2dual)
208 {
209 ENUMERATE_ITEM(item, family->allItems().own())
210 {
211 auto connected_items = accessor.connectedItems(ItemLocalId(item)) ;
212 if(connected_items.size()>0)
213 {
214 dofs_new_owner[*item] = dualitems_new_owner[connected_items[0]] ;
215 }
216 }
217 }
218 else
219 {
220 ENUMERATE_ITEM(item, dualitem_family->allItems())
221 {
222 ENUMERATE_ITEM(connected_item,accessor.connectedItems(ItemLocalId(item)))
223 {
224 dofs_new_owner[*connected_item] = dualitems_new_owner[*item] ;
225 }
226 }
227 }
228 }
229 }
230 }
231 }
232 // Dof with DoF connectivity
233 IItemFamily* links_family = m_mesh->findItemFamily(IK_DoF, mesh::GraphDoFs::linkFamilyName(), false);
234 if(links_family)
235 {
236 VariableItemInt32& links_new_owner(links_family->itemsNewOwner());
237 IItemFamily* dualnodes_family = m_mesh->findItemFamily(IK_DoF, mesh::GraphDoFs::dualNodeFamilyName(), false);
238 if(dualnodes_family)
239 {
240 VariableItemInt32& dualnodes_new_owner(dualnodes_family->itemsNewOwner());
241 auto connectivity_name = mesh::connectivityName(links_family,dualnodes_family) ;
242 auto connectivity = m_mesh->itemFamilyNetwork()->getConnectivity(links_family,dualnodes_family,connectivity_name) ;
243 if(connectivity)
244 {
245 ConnectivityItemVector accessor(connectivity);
246 ENUMERATE_ITEM(item, links_family->allItems().own())
247 {
248 auto connected_items = accessor.connectedItems(ItemLocalId(item)) ;
249 if(connected_items.size()>0)
250 {
251 links_new_owner[*item] = dualnodes_new_owner[connected_items[0]] ;
252 }
253 }
254 }
255 }
256 }
257 }
258}
259
260/*---------------------------------------------------------------------------*/
261/*---------------------------------------------------------------------------*/
262
263void UnstructuredMeshUtilities::
264localIdsFromConnectivity(eItemKind item_kind,
265 IntegerConstArrayView items_nb_node,
266 Int64ConstArrayView items_connectivity,
267 Int32ArrayView local_ids,
268 bool allow_null)
269{
270 Integer nb_item = items_nb_node.size();
271 if (item_kind!=IK_Face && item_kind!=IK_Cell)
272 throw InvalidArgumentException(A_FUNCINFO,"item_kind",
273 "IK_Cell or IK_Face expected",
274 (int)item_kind);
275 if (item_kind==IK_Cell)
276 throw NotImplementedException(A_FUNCINFO,"not implemented for cells");
277
278 if (nb_item!=local_ids.size())
279 throw InvalidArgumentException(A_FUNCINFO,"local_ids",
280 "Size different from 'items_nb_node'",
281 (int)item_kind);
282
283 Integer item_connectivity_index = 0;
285 buf.reserve(256);
286 ItemInternalList nodes(m_mesh->itemsInternal(IK_Node));
287 for( Integer i=0; i<nb_item; ++i ){
288 Integer current_nb_node = items_nb_node[i];
289 Int64ConstArrayView current_nodes(current_nb_node,items_connectivity.data()+item_connectivity_index);
290 item_connectivity_index += current_nb_node;
291 if (item_kind==IK_Face){
292 buf.resize(current_nb_node);
293 mesh_utils::reorderNodesOfFace(current_nodes,buf);
294 Int64 first_node_uid = buf[0];
295 Int64ArrayView first_node_uid_array(1,&first_node_uid);
296 Int32 first_node_lid = CheckedConvert::toInt32(buf[0]);
297 Int32ArrayView first_node_lid_array(1,&first_node_lid);
298 m_mesh->nodeFamily()->itemsUniqueIdToLocalId(first_node_lid_array,first_node_uid_array,!allow_null);
299 if (first_node_lid == NULL_ITEM_LOCAL_ID){
300 if (allow_null){
301 local_ids[i] = NULL_ITEM_LOCAL_ID;
302 }
303 else {
304 StringBuilder sb("Face with nodes (");
305 for(Integer j=0;j<current_nb_node;++j) {
306 if (j != 0) sb += " ";
307 sb += current_nodes[j];
308 }
309 sb += ") not found (first node ";
310 sb += first_node_uid;
311 sb += " not found)";
313 }
314 continue;
315 }
316
317 Node node(nodes[first_node_lid]);
318 Face face(mesh_utils::getFaceFromNodesUnique(node,buf));
319 if (face.null()){
320 if (allow_null){
321 local_ids[i] = NULL_ITEM_LOCAL_ID;
322 }
323 else {
324 StringBuilder sb("Face with nodes (");
325 for(Integer j=0;j<current_nb_node;++j) {
326 if (j != 0) sb += " ";
327 sb += current_nodes[j];
328 }
329 sb += ") not found";
331 }
332 }
333 else
334 local_ids[i] = face.localId();
335 }
336 }
337}
338
339/*---------------------------------------------------------------------------*/
340/*---------------------------------------------------------------------------*/
341
342Real3 UnstructuredMeshUtilities::
343_round(Real3 value)
344{
345 Real3 rvalue = value;
346 // Evite les problèmes d'arrondi numérique
347 if (math::isNearlyZero(value.x))
348 rvalue.x = 0.0;
349 if (math::isNearlyZero(value.y))
350 rvalue.y = 0.0;
351 if (math::isNearlyZero(value.z))
352 rvalue.z = 0.0;
353
354 if (math::isNearlyEqual(value.x,1.0))
355 rvalue.x = 1.0;
356 if (math::isNearlyEqual(value.y,1.0))
357 rvalue.y = 1.0;
358 if (math::isNearlyEqual(value.z,1.0))
359 rvalue.z = 1.0;
360
361 if (math::isNearlyEqual(value.x,-1.0))
362 rvalue.x = -1.0;
363 if (math::isNearlyEqual(value.y,-1.0))
364 rvalue.y = -1.0;
365 if (math::isNearlyEqual(value.z,-1.0))
366 rvalue.z = -1.0;
367
368 return rvalue;
369}
370
371/*---------------------------------------------------------------------------*/
372/*---------------------------------------------------------------------------*/
373
374Real3 UnstructuredMeshUtilities::
375computeNormal(const FaceGroup& face_group,const VariableNodeReal3& nodes_coord)
376{
377 String func_name = "UnstructuredMeshUtilities::computeNormal";
378 IParallelMng* pm = m_mesh->parallelMng();
379 Integer rank = pm->commRank();
380 //Integer nb_item = face_group.size();
381 Integer nb_node = 0;
382 Integer mesh_dimension = m_mesh->dimension();
383 // 'global_normal' et 'total_global_normal' servent à determiner
384 // une direction pour la surface en supposant qu'il s'agit
385 // d'une surface externe. La direction sert ensuite à orienter
386 // la normale calculée dans le sens normal sortante.
387 Real3 global_normal = Real3::null();
388 Real nb_global_normal_used_node = 0.0;
389 ENUMERATE_FACE(iface,face_group){
390 Face face = *iface;
391 Integer face_nb_node = face.nbNode();
392 nb_node += face_nb_node;
393
394 if (face.isSubDomainBoundary() && face.isOwn()){
395 Real3 face_normal;
396 if (mesh_dimension==3){
397 Real3 v1 = nodes_coord[face.node(1)] - nodes_coord[face.node(0)];
398 Real3 v2 = nodes_coord[face.node(2)] - nodes_coord[face.node(0)];
399 face_normal = math::vecMul(v1,v2);
400 }
401 else if (mesh_dimension==2){
402 Real3 dir = nodes_coord[face.node(0)] - nodes_coord[face.node(1)];
403 face_normal.x = - dir.y;
404 face_normal.y = dir.x;
405 face_normal.z = 0.0;
406 }
407 if (face.boundaryCell()==face.frontCell())
408 face_normal = -face_normal;
409 //info() << "ADD NORMAL normal=" << face_normal;
410 global_normal = global_normal + face_normal;
411 nb_global_normal_used_node += 1.0;
412 }
413 }
414 // L'algorithme tente de déterminer les noeuds aux extrémités de la surface
415 // On considère qu'il s'agit de ceux qui n'appartiennent qu'à une seule
416 // face de la surface. Cela fonctionne bien si toutes les faces ont au moins
417 // quatres arêtes. S'il y a des triangles, on suppose que ce n'est pas
418 // partout et dans ce cas on a au moins deux noeuds extrèmes connus.
419 typedef HashTableMapT<Int32,Int32> NodeOccurenceMap;
420 // Range dans la table pour chaque noeud le nombre de faces de la surface
421 // auxquelles il est connecté.
422 NodeOccurenceMap nodes_occurence(nb_node*2,true,nb_node);
423 ENUMERATE_FACE(iface,face_group){
424 Face face = *iface;
425 for( NodeLocalId inode : face.nodeIds() ){
426 NodeOccurenceMap::Data* v = nodes_occurence.lookup(inode);
427 if (v)
428 ++v->value();
429 else{
430 //info() << " ADD NODE " << ItemPrinter(*inode) << " coord=" << nodes_coord[inode];
431 nodes_occurence.add(inode,1);
432 }
433 }
434 }
435 UniqueArray<Node> single_nodes;
436 NodeLocalIdToNodeConverter nodes_internal(m_mesh->nodeFamily());
437 nodes_occurence.each([&](NodeOccurenceMap::Data* d){
438 if (d->value()==1){
439 Node node = nodes_internal[d->key()];
440 // En parallèle, ne traite que les noeuds qui nous appartiennent
441 if (node.owner()==rank){
442 single_nodes.add(node);
443 //info() << "SINGLE NODE OWNER lid=" << d->key() << " " << ItemPrinter(node)
444 // << " coord=" << nodes_coord[node];
445 }
446 }
447 });
448 // Chaque sous-domaine collecte les coordonnées des noeuds des autres sous-domaines
449 Integer nb_single_node = single_nodes.size();
450 //info() << "NB SINGLE NODE= " << nb_single_node;
451 Integer total_nb_single_node = pm->reduce(Parallel::ReduceSum,nb_single_node);
452 Real3 total_global_normal;
453 {
454 Real all_min = 0.0;
455 Real all_max = 0.0;
456 Real all_sum = 0.0;
457 Int32 min_rank = -1;
458 Int32 max_rank = -1;
459 pm->computeMinMaxSum(nb_global_normal_used_node,all_min,all_max,all_sum,min_rank,max_rank);
460 Real3 buf;
461 if (max_rank==rank){
462 // Je suis celui qui a le point le plus loin. Je l'envoie aux autres.
463 buf = global_normal;
464 }
465 pm->broadcast(Real3ArrayView(1,&buf),max_rank);
466 total_global_normal = buf;
467 }
468
469 info() << "TOTAL SINGLE NODE= " << total_nb_single_node << " surface=" << face_group.name()
470 << " global_normal = " << total_global_normal;
471 if (total_nb_single_node<2){
472 ARCANE_FATAL("not enough nodes connected to only one face");
473 }
474 Real3UniqueArray coords(nb_single_node);
475 for( Integer i=0; i<nb_single_node; ++i ){
476 coords[i] = nodes_coord[single_nodes[i]];
477 }
478 //Array<Real> all_nodes_coord_real;
479 UniqueArray<Real3> all_nodes_coord;
480 pm->allGatherVariable(coords,all_nodes_coord);
481 //info() << " ALL NODES COORD n=" << all_nodes_coord_real.size();
482 //Array<Real3> all_nodes_coord(total_nb_single_node);
483 Real3 barycentre = Real3::null();
484 for( Integer i=0; i<total_nb_single_node; ++i ){
485 barycentre += all_nodes_coord[i];
486 }
487 barycentre /= (Real)total_nb_single_node;
488 if (total_nb_single_node==2 && m_mesh->dimension()==3){
489 // On a que deux noeuds. Il en faut au moins un troisième pour déterminer
490 // la normale au plan. Pour cela, chaque processeur cherche le noeud de
491 // la surface qui est le plus éloigné du barycentre des deux noeuds déjà
492 // trouvé. Ce noeud le plus éloigné servira de troisième point pour le
493 // plan.
494 Real max_distance = 0.0;
495 Node farthest_node;
496 Real3 s0 = all_nodes_coord[0];
497 Real3 s1 = all_nodes_coord[1];
498 nodes_occurence.each([&](NodeOccurenceMap::Data* d){
499 Node node = nodes_internal[d->key()];
500 Real3 coord = nodes_coord[node];
501 // Ne traite pas les deux noeuds du déjà trouvés
502 if (math::isNearlyEqual(coord,s0))
503 return;
504 if (math::isNearlyEqual(coord,s1))
505 return;
506 Real distance = (coord - barycentre).squareNormL2();
507 if (distance>max_distance){
508 Real3 normal = math::cross(coord-s0, s1-s0);
509 // On ne prend le noeud que s'il n'est pas aligné avec les deux autres.
510 if (!math::isNearlyZero(normal.squareNormL2())){
511 max_distance = distance;
512 farthest_node = node;
513 }
514 }
515 });
516
517 if (!farthest_node.null()){
518 info() << " FARTHEST NODE= " << ItemPrinter(farthest_node) << " dist=" << max_distance;
519 }
520 {
521 Real3 farthest_coord = _broadcastFarthestNode(max_distance,farthest_node,nodes_coord);
522 info() << " FARTHEST NODE ALL coord=" << farthest_coord;
523 all_nodes_coord.add(farthest_coord);
524 }
525 }
526 // Trie les noeuds pour que le calcul ne dépende pas de l'ordre des
527 // opérations en parallèle.
528 std::sort(std::begin(all_nodes_coord),std::end(all_nodes_coord));
529 Integer nb_final_node = all_nodes_coord.size();
530 Real3 full_normal = Real3::null();
531 info() << " NB FINAL NODE=" << nb_final_node;
532 for( Integer i=0; i<nb_final_node; ++i )
533 info() << " NODE=" << all_nodes_coord[i];
534 if (m_mesh->dimension()==2){
535 if (nb_final_node!=2){
536 ARCANE_FATAL("should have 2 border nodes in 2D mesh");
537 }
538 Real3 direction = all_nodes_coord[1] - all_nodes_coord[0];
539 full_normal.x = - direction.y;
540 full_normal.y = direction.x;
541 }
542 else if (m_mesh->dimension()==3){
543 //nb_final_node = 3; // On prend que les 3 premiers points.
544 // NOTE: on pourrait prendre tous les points, car si les trois premiers sont
545 // alignés, la normale ne sera pas bonne.
546 // Si on prend tous les points, il faut être sur qu'ils soient ordonnés
547 // de telle sorte que la normale de trois points consécutifs est toujours
548 // dans le même sens. Cela signifie si on a quatres points par exemple,
549 // que ces 4 points forment un quadrangle non croisé (pas en forme
550 // de papillon)
551 Real3 first_normal = math::vecMul(all_nodes_coord[2]-all_nodes_coord[0],
552 all_nodes_coord[1]-all_nodes_coord[0]);
553 for( Integer i=0; i<nb_final_node; ++i ){
554 Real3 s0 = all_nodes_coord[i];
555 Real3 s1 = all_nodes_coord[(i+nb_final_node-1)%nb_final_node];
556 Real3 s2 = all_nodes_coord[(i+1)%nb_final_node];
557 Real3 normal = math::vecMul(s2-s0, s1-s0);
558 info() << " ADD NORMAL: " << normal;
559 full_normal += normal;
560 info() << " FULL: " << full_normal;
561 }
562 if (math::isNearlyZero(full_normal.squareNormL2()))
563 full_normal = first_normal;
564 }
565 else
566 ARCANE_FATAL("invalid mesh dimension (should be 2 or 3)");
567 Real a = full_normal.normL2();
568 if (math::isZero(a))
569 ARCANE_FATAL("invalid value for normal");
570 full_normal /= a;
571 Real b = total_global_normal.normL2();
572 Real dir = 0.0;
573 info() << " Normal is " << full_normal;
574 if (!math::isZero(b)){
575 total_global_normal /= b;
576 dir = math::scaMul(full_normal,total_global_normal);
577 if (dir<0.0)
578 full_normal = -full_normal;
579 }
580 full_normal = _round(full_normal);
581 info() << " Final normal is " << full_normal << " dir=" << dir;
582 return full_normal;
583}
584
585/*---------------------------------------------------------------------------*/
586/*---------------------------------------------------------------------------*/
594Real3 UnstructuredMeshUtilities::
595computeDirection(const NodeGroup& node_group,
596 const VariableNodeReal3& nodes_coord,
597 Real3* n1,Real3* n2)
598{
599 String func_name = "UnstructuredMeshUtilities::computeDirection";
600 IParallelMng* pm = m_mesh->parallelMng();
601 //Integer rank = pm->commRank();
602 //Integer nb_node_own = node_group.own().size();
603 //Integer total_nb_node_own = pm->reduce(Parallel::ReduceSum,nb_node_own);
604 Real3 barycentre = Real3::null();
605 ENUMERATE_NODE(inode,node_group.own()){
606 barycentre += nodes_coord[inode];
607 }
608 Real r_barycentre[3];
609 r_barycentre[0] = barycentre.x;
610 r_barycentre[1] = barycentre.y;
611 r_barycentre[2] = barycentre.z;
612 RealArrayView rav(3,r_barycentre);
614 barycentre = Real3(rav[0],rav[1],rav[2]);
615 debug() << " BARYCENTRE COMPUTE DIRECTION = " << barycentre;
616 pm->barrier();
617
618 // Détermine le noeud le plus éloigné du barycentre
619 Real3 first_boundary_coord;
620 {
621 Real max_distance = 0.0;
622 Node farthest_node;
623 ENUMERATE_NODE(inode,node_group.own()){
624 Node node = *inode;
625 Real3 coord = nodes_coord[node];
626 Real distance = (coord - barycentre).squareNormL2();
627 if (distance>max_distance){
628 max_distance = distance;
629 farthest_node = node;
630 }
631 }
632 first_boundary_coord = _broadcastFarthestNode(max_distance,farthest_node,nodes_coord);
633 if (n1)
634 *n1 = first_boundary_coord;
635 }
636 Real3 second_boundary_coord;
637 {
638 Real max_distance = 0.0;
639 Node farthest_node;
640 ENUMERATE_NODE(inode,node_group.own()){
641 Node node = *inode;
642 Real3 coord = nodes_coord[node];
643 Real distance = (coord - first_boundary_coord).squareNormL2();
644 if (distance>max_distance){
645 max_distance = distance;
646 farthest_node = node;
647 }
648 }
649 second_boundary_coord = _broadcastFarthestNode(max_distance,farthest_node,nodes_coord);
650 if (n2)
651 *n2 = second_boundary_coord;
652 }
653 Real3 direction = second_boundary_coord - first_boundary_coord;
654 Real norm = direction.normL2();
655 if (math::isZero(norm)){
656 ARCANE_FATAL("Direction is null for group '{0}' first_coord={1} second_coord={2}",
657 node_group.name(),first_boundary_coord,second_boundary_coord);
658 }
659 return direction / norm;
660}
661
662/*---------------------------------------------------------------------------*/
663/*---------------------------------------------------------------------------*/
664
665Real3 UnstructuredMeshUtilities::
666_broadcastFarthestNode(Real distance,const Node& farthest_node,
667 const VariableNodeReal3& nodes_coord)
668{
669 String func_name = "UnstructuredMeshUtilities::_broadcastFarthestNode()";
670 IParallelMng* pm = m_mesh->parallelMng();
671 Integer rank = pm->commRank();
672
673 Real all_min_distance = 0.0;
674 Real all_max_distance = 0.0;
675 Real all_sum_distance = 0.0;
676 Int32 min_rank = -1;
677 Int32 max_rank = -1;
678 pm->computeMinMaxSum(distance,all_min_distance,all_max_distance,all_sum_distance,min_rank,max_rank);
679 Real3 buf;
680 if (!farthest_node.null())
681 debug() << " FARTHEST NODE myself coord=" << nodes_coord[farthest_node]
682 << " distance=" << distance
683 << " rank=" << max_rank
684 << " max_distance=" << all_max_distance;
685 else
686 debug() << " FARTHEST NODE myself coord=none"
687 << " distance=" << distance
688 << " rank=" << max_rank
689 << " max_distance=" << all_max_distance;
690
691 if (max_rank==rank){
692 if (farthest_node.null())
693 ARCANE_FATAL("can not find farthest node");
694 // Je suis celui qui a le point le plus loin. Je l'envoie aux autres.
695 buf = nodes_coord[farthest_node];
696 debug() << " I AM FARTHEST NODE ALL coord=" << buf;
697 }
698 pm->broadcast(Real3ArrayView(1,&buf),max_rank);
699 Real3 farthest_coord(buf);
700 debug() << " FARTHEST NODE ALL coord=" << farthest_coord
701 << " rank=" << max_rank
702 << " max_distance=" << all_max_distance;
703 return farthest_coord;
704}
705
706/*---------------------------------------------------------------------------*/
707/*---------------------------------------------------------------------------*/
708
709void UnstructuredMeshUtilities::
710computeAdjency(ItemPairGroup adjency_array,eItemKind link_kind,Integer nb_layer)
711{
712 m_compute_adjacency_functor->computeAdjacency(adjency_array, link_kind, nb_layer);
713}
714
715/*---------------------------------------------------------------------------*/
716/*---------------------------------------------------------------------------*/
717
718bool UnstructuredMeshUtilities::
719writeToFile(const String& file_name,const String& service_name)
720{
721 ServiceBuilder<IMeshWriter> sb(m_mesh->handle());
722 auto mesh_writer = sb.createReference(service_name,SB_AllowNull);
723
724 if (!mesh_writer){
725 UniqueArray<String> available_names;
726 sb.getServicesNames(available_names);
727 warning() << String::format("The specified service '{0}' to write the mesh is not available."
728 " Valid names are {1}",service_name,available_names);
729 return true;
730 }
731 mesh_writer->writeMeshToFile(m_mesh,file_name);
732 return false;
733}
734
735/*---------------------------------------------------------------------------*/
736/*---------------------------------------------------------------------------*/
737
738void UnstructuredMeshUtilities::
739partitionAndExchangeMeshWithReplication(IMeshPartitionerBase* partitioner,
740 bool initial_partition)
741{
742 IPrimaryMesh* primary_mesh = partitioner->primaryMesh();
743 IMesh* mesh = primary_mesh;
744 if (mesh!=this->m_mesh)
745 throw ArgumentException(A_FUNCINFO,"partitioner->mesh() != this->m_mesh");
746
747 IParallelMng* pm = mesh->parallelMng();
748 ITimeStats* ts = pm->timeStats();
750 bool has_replication = pr->hasReplication();
751
752 // En mode réplication, seul le premier réplicat fait l'équilibrage.
753 // Il doit ensuite envoyer aux autres ces informations de maillage
754 // pour qu'ils aient le même maillage.
755
756 info() << "Partition start date=" << platform::getCurrentDateTime();
757 if (pr->isMasterRank()){
758 Timer::Action ts_action1(ts,"MeshPartition",true);
759 info() << "Partitioning the mesh (initial?=" << initial_partition << ")";
760 partitioner->partitionMesh(initial_partition);
761 }
762 else
763 info() << "Waiting for partition information from the master replica";
764
765 if (has_replication){
766 pm->barrier();
767
768 // Vérifie que toute les familles sont les mêmes.
769 mesh->checker()->checkValidReplication();
770
771 Int32 replica_master_rank = pr->masterReplicationRank();
772 IParallelMng* rep_pm = pr->replicaParallelMng();
773
774 // Seul le replica maitre a les bons propriétaires. Il faut mettre
775 // à jour les autres à partir de celui-ci. Pour cela on synchronize
776 // les propriétaires des mailles et ensuite on met à jour les autres familles
777 // à partir des propriétaires des mailles.
778 {
779 Int32ArrayView owners = mesh->cellFamily()->itemsNewOwner().asArray();
780 rep_pm->broadcast(owners,replica_master_rank);
781 if (!pr->isMasterRank()){
783 }
784 }
785 }
786 info() << "Partition end date=" << platform::getCurrentDateTime();
787 {
788 Timer::Action ts_action2(ts,"MeshExchange",true);
789 primary_mesh->exchangeItems();
790 }
791 info() << "Exchange end date=" << platform::getCurrentDateTime();
792 partitioner->notifyEndPartition();
793
794 // Il faut recompacter pour se retrouver dans la même situation que
795 // si on avait juste lu un maillage directement (qui fait un prepareForDump()
796 // lors de l'appel à endAllocate()).
797 // On le fait aussi en cas de réplication pour éviter d'éventuelles
798 // incohérences si par la suite certains réplicas appellent cette
799 // méthode et pas les autres.
800 // TODO: regarder s'il faut le faire aussi en cas de repartitionnement.
801 if (initial_partition || has_replication)
802 mesh->prepareForDump();
803}
804
805/*---------------------------------------------------------------------------*/
806/*---------------------------------------------------------------------------*/
807
808void UnstructuredMeshUtilities::
809mergeNodes(Int32ConstArrayView nodes_local_id,
810 Int32ConstArrayView nodes_to_merge_local_id)
811{
812 mesh::MeshNodeMerger merger(m_mesh);
813 merger.mergeNodes(nodes_local_id,nodes_to_merge_local_id);
814}
815
816/*---------------------------------------------------------------------------*/
817/*---------------------------------------------------------------------------*/
818
819void UnstructuredMeshUtilities::
820computeAndSetOwnersForNodes()
821{
822 mesh::ItemsOwnerBuilder owner_builder(m_mesh);
823 owner_builder.computeNodesOwner();
824}
825
826/*---------------------------------------------------------------------------*/
827/*---------------------------------------------------------------------------*/
828
829void UnstructuredMeshUtilities::
830computeAndSetOwnersForFaces()
831{
832 mesh::ItemsOwnerBuilder owner_builder(m_mesh);
833 owner_builder.computeFacesOwner();
834}
835
836
837/*---------------------------------------------------------------------------*/
838/*---------------------------------------------------------------------------*/
839namespace
840{
841 void _recomputeUniqueIds(IItemFamily* family)
842 {
843 ITraceMng* tm = family->traceMng();
844 eItemKind ik = family->itemKind();
845 SmallArray<Int64> unique_ids;
846 ENUMERATE_ (ItemWithNodes, iitem, family->allItems()) {
847 ItemWithNodes item = *iitem;
848 Int32 index = 0;
849 unique_ids.resize(item.nbNode());
850 for (Node node : item.nodes()) {
851 unique_ids[index] = node.uniqueId();
852 ++index;
853 }
854 Int64 new_uid = MeshUtils::generateHashUniqueId(unique_ids);
855 //if (ik==IK_Face)
856 //tm->info() << "Face uid=" << item.uniqueId() << " nb_node=" << item.nbNode()
857 // << " node0=" << item.node(0).uniqueId() << " node1=" << item.node(1).uniqueId();
858 item.mutableItemBase().setUniqueId(new_uid);
859 }
861 }
862} // namespace
863
864/*---------------------------------------------------------------------------*/
865/*---------------------------------------------------------------------------*/
866
867void UnstructuredMeshUtilities::
868recomputeItemsUniqueIdFromNodesUniqueId()
869{
870 IMesh* mesh = m_mesh;
872 ITraceMng* tm = mesh->traceMng();
873
874 tm->info() << "Calling RecomputeItemsUniqueIdFromNodesUniqueId()";
875 // D'abord indiquer que les noeuds ont changés pour éventuellement
876 // remettre à jour l'orientation des faces.
877 mesh->nodeFamily()->notifyItemsUniqueIdChanged();
878 _recomputeUniqueIds(mesh->edgeFamily());
879 _recomputeUniqueIds(mesh->faceFamily());
880 _recomputeUniqueIds(mesh->cellFamily());
881}
882
883/*---------------------------------------------------------------------------*/
884/*---------------------------------------------------------------------------*/
885
886} // End namespace Arcane
887
888/*---------------------------------------------------------------------------*/
889/*---------------------------------------------------------------------------*/
#define ARCANE_CHECK_POINTER(ptr)
Macro retournant le pointeur ptr s'il est non nul ou lancant une exception s'il est nul.
#define ARCANE_FATAL(...)
Macro envoyant une exception FatalErrorException.
#define ENUMERATE_FACE(name, group)
Enumérateur générique d'un groupe de faces.
#define ENUMERATE_(type, name, group)
Enumérateur générique d'un groupe d'entité
#define ENUMERATE_PARTICLE(name, group)
Enumérateur générique d'un groupe de particules.
#define ENUMERATE_EDGE(name, group)
Enumérateur générique d'un groupe d'arêtes.
#define ENUMERATE_ITEM(name, group)
Enumérateur générique d'un groupe de noeuds.
#define ENUMERATE_NODE(name, group)
Enumérateur générique d'un groupe de noeuds.
bool reorderNodesOfFace(Int64ConstArrayView before_ids, Int64ArrayView after_ids)
Réordonne les noeuds d'une face.
Integer size() const
Nombre d'éléments du vecteur.
Exception lorsqu'un argument est invalide.
constexpr Integer size() const noexcept
Retourne la taille du tableau.
void resize(Int64 s)
Change le nombre d'éléments du tableau à s.
void reserve(Int64 new_capacity)
Réserve le mémoire pour new_capacity éléments.
void add(ConstReferenceType val)
Ajoute l'élément val à la fin du tableau.
Maille d'un maillage.
Definition Item.h:1191
Gère la récupération des informations de connectivité.
ItemVectorView connectedItems(ItemLocalId item)
Retourne les entités connectées à item.
constexpr const_pointer data() const noexcept
Pointeur sur la mémoire allouée.
constexpr Integer size() const noexcept
Nombre d'éléments du tableau.
Arête d'une maille.
Definition Item.h:809
Face d'une maille.
Definition Item.h:944
Cell frontCell() const
Maille devant la face (maille nulle si aucune)
Definition Item.h:1620
bool isSubDomainBoundary() const
Indique si la face est au bord du sous-domaine (i.e nbCell()==1)
Definition Item.h:1038
Cell boundaryCell() const
Maille associée à cette face frontière (maille nulle si aucune)
Definition Item.h:1608
Table de hachage pour tableaux associatifs.
Interface d'une famille d'entités.
Definition IItemFamily.h:84
virtual void notifyItemsUniqueIdChanged()=0
Notifie que les numéros uniques des entités ont été modifiées.
virtual ItemGroup allItems() const =0
Groupe de toutes les entités.
virtual eItemKind itemKind() const =0
Genre des entités.
virtual ITraceMng * traceMng() const =0
Gestionnaire de trace associé
virtual VariableItemInt32 & itemsNewOwner()=0
Variable contenant le numéro du nouveau sous-domaine propriétaire de l'entité.
Interface d'un partitionneur de maillage.
virtual IPrimaryMesh * primaryMesh()=0
Maillage associé
virtual void partitionMesh(bool initial_partition)=0
virtual void notifyEndPartition()=0
Notification lors de la fin d'un repartionnement (après échange des entités)
Interface du gestionnaire de parallélisme pour un sous-domaine.
virtual void computeMinMaxSum(char val, char &min_val, char &max_val, char &sum_val, Int32 &min_rank, Int32 &max_rank)=0
Calcule en une opération la somme, le min, le max d'une valeur.
virtual Int32 commRank() const =0
Rang de cette instance dans le communicateur.
virtual IParallelReplication * replication() const =0
Informations sur la réplication.
virtual ITimeStats * timeStats() const =0
Gestionnaire de statistiques associé (peut être nul)
virtual void allGatherVariable(ConstArrayView< char > send_buf, Array< char > &recv_buf)=0
Effectue un regroupement sur tous les processeurs.
virtual char reduce(eReduceType rt, char v)=0
Effectue la réduction de type rt sur le réel v et retourne la valeur.
virtual void barrier()=0
Effectue une barière.
Informations sur la réplication des sous-domaines en parallèle.
virtual bool hasReplication() const =0
Indique si la réplication est active.
virtual IParallelMng * replicaParallelMng() const =0
Communicateur associé à tous les réplicats représentant un même sous-domaine.
virtual Int32 masterReplicationRank() const =0
Rang dans la réplication du maître.
virtual bool isMasterRank() const =0
Indique si ce rang de réplication est le maître.
virtual void exchangeItems()=0
Change les sous-domaines propriétaires des entités.
Interface gérant les statistiques sur les temps d'exécution.
Definition ITimeStats.h:43
Interface du gestionnaire de traces.
virtual TraceMessage info()=0
Flot pour un message d'information.
Exception lorsqu'une erreur fatale est survenue.
const String & name() const
Nom du groupe.
Definition ItemGroup.h:76
ItemGroup own() const
Groupe équivalent à celui-ci mais contenant uniquement les éléments propres au sous-domaine.
Definition ItemGroup.cc:189
Index d'un Item dans une variable.
Definition ItemLocalId.h:41
Tableau de listes d'entités.
Classe utilitaire pour imprimer les infos sur une entité.
Definition ItemPrinter.h:35
Elément de maillage s'appuyant sur des noeuds (Edge,Face,Cell).
Definition Item.h:724
Node node(Int32 i) const
i-ème noeud de l'entité
Definition Item.h:779
NodeConnectedListViewType nodes() const
Liste des noeuds de l'entité
Definition Item.h:782
Int32 nbNode() const
Nombre de noeuds de l'entité
Definition Item.h:776
NodeLocalIdView nodeIds() const
Liste des noeuds de l'entité
Definition Item.h:785
impl::MutableItemBase mutableItemBase() const
Partie interne modifiable de l'entité.
Definition Item.h:374
constexpr Int32 localId() const
Identifiant local de l'entité dans le sous-domaine du processeur.
Definition Item.h:219
bool isOwn() const
true si l'entité est appartient au sous-domaine
Definition Item.h:253
Int32 owner() const
Numéro du sous-domaine propriétaire de l'entité
Definition Item.h:238
constexpr bool null() const
true si l'entité est nul (i.e. non connecté au maillage)
Definition Item.h:216
Classe pour convertir un NodeLocalId vers une arête.
Noeud d'un maillage.
Definition Item.h:573
Exception lorsqu'une fonction n'est pas implémentée.
Particule.
Definition Item.h:1397
Cell cell() const
Maille à laquelle appartient la particule. Il faut appeler setCell() avant d'appeler cette fonction....
Definition Item.h:1459
Classe gérant un vecteur de réel de dimension 3.
Definition Real3.h:132
__host__ __device__ Real normL2() const
Retourne la norme L2 du triplet .
Definition Real3.h:521
constexpr __host__ __device__ Real squareNormL2() const
Retourne la norme L2 au carré du triplet .
Definition Real3.h:408
Classe utilitaire pour instantier un service d'une interface donnée.
Ref< InterfaceType > createReference(const String &name, eServiceBuilderProperties properties=SB_None)
Créé une instance implémentant l'interface InterfaceType.
void getServicesNames(Array< String > &names) const
Remplit names avec les noms des services disponibles pour instantier cette interface.
Tableau 1D de données avec buffer pré-alloué sur la pile.
Definition SmallArray.h:89
Constructeur de chaîne de caractère unicode.
String toString() const
Retourne la chaîne de caractères construite.
Chaîne de caractères unicode.
Postionne le nom de l'action en cours d'exécution.
Definition Timer.h:110
TraceMessageDbg debug(Trace::eDebugLevel=Trace::Medium) const
Flot pour un message de debug.
TraceMessage info() const
Flot pour un message d'information.
TraceMessage warning() const
Flot pour un message d'avertissement.
Vecteur 1D de données avec sémantique par valeur (style STL).
void changeOwnersFromCells() override
Positionne les nouveaux propriétaires des noeuds, arêtes et faces à partir des mailles.
Classe générique pour calculer les propriétaires des entités.
__host__ __device__ Real3 vecMul(Real3 u, Real3 v)
Produit vectoriel de u par v. dans .
Definition MathUtils.h:52
__host__ __device__ Real scaMul(Real2 u, Real2 v)
Produit scalaire de u par v dans .
Definition MathUtils.h:113
__host__ __device__ Real3 cross(Real3 v1, Real3 v2)
Produit vectoriel de deux vecteurs à 3 composantes.
Definition MathUtils.h:723
ItemGroupT< Face > FaceGroup
Groupe de faces.
Definition ItemTypes.h:178
ItemGroupT< Node > NodeGroup
Groupe de noeuds.
Definition ItemTypes.h:167
MeshVariableScalarRefT< Node, Real3 > VariableNodeReal3
Grandeur au noeud de type coordonnées.
ItemVariableScalarRefT< Int32 > VariableItemInt32
Grandeur de type entier 32 bits.
Int32 toInt32(Int64 v)
Converti un Int64 en un Int32.
@ ReduceSum
Somme des valeurs.
constexpr __host__ __device__ bool isNearlyEqual(const _Type &a, const _Type &b)
Teste si deux valeurs sont à un peu près égales. Pour les types entiers, cette fonction est équivalen...
Definition Numeric.h:212
bool isZero(const BuiltInProxy< _Type > &a)
Teste si une valeur est exactement égale à zéro.
ARCCORE_BASE_EXPORT String getCurrentDateTime()
Date et l'heure courante sous la forme ISO 8601.
-*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
ArrayView< Int64 > Int64ArrayView
Equivalent C d'un tableau à une dimension d'entiers 64 bits.
Definition UtilsTypes.h:538
@ SB_AllowNull
Autorise l'absence du service.
UniqueArray< Int64 > Int64UniqueArray
Tableau dynamique à une dimension d'entiers 64 bits.
Definition UtilsTypes.h:426
std::int64_t Int64
Type entier signé sur 64 bits.
ArrayView< Real3 > Real3ArrayView
Equivalent C d'un tableau à une dimension de Real3.
Definition UtilsTypes.h:554
Int32 Integer
Type représentant un entier.
UniqueArray< Real3 > Real3UniqueArray
Tableau dynamique à une dimension de vecteurs de rang 3.
Definition UtilsTypes.h:450
ConstArrayView< Int32 > Int32ConstArrayView
Equivalent C d'un tableau à une dimension d'entiers 32 bits.
Definition UtilsTypes.h:569
ConstArrayView< ItemInternal * > ItemInternalList
Type de la liste interne des entités.
Definition ItemTypes.h:466
ConstArrayView< Int64 > Int64ConstArrayView
Equivalent C d'un tableau à une dimension d'entiers 64 bits.
Definition UtilsTypes.h:567
ArrayView< Int32 > Int32ArrayView
Equivalent C d'un tableau à une dimension d'entiers 32 bits.
Definition UtilsTypes.h:540
eItemKind
Genre d'entité de maillage.
@ IK_Particle
Entité de maillage de genre particule.
@ IK_Node
Entité de maillage de genre noeud.
@ IK_Cell
Entité de maillage de genre maille.
@ IK_Face
Entité de maillage de genre face.
@ IK_DoF
Entité de maillage de genre degre de liberte.
@ IK_Edge
Entité de maillage de genre arête.
double Real
Type représentant un réel.
ConstArrayView< Integer > IntegerConstArrayView
Equivalent C d'un tableau à une dimension d'entiers.
Definition UtilsTypes.h:573
ArrayView< Real > RealArrayView
Equivalent C d'un tableau à une dimension de réels.
Definition UtilsTypes.h:546
std::int32_t Int32
Type entier signé sur 32 bits.
Real y
deuxième composante du triplet
Definition Real3.h:36
Real z
troisième composante du triplet
Definition Real3.h:37
Real x
première composante du triplet
Definition Real3.h:35