OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
subdivision.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 #include "subdivision.h"
36 #include "material.h"
37 #include "errorestimator.h"
38 #include "domain.h"
39 #include "node.h"
40 #include "element.h"
41 #include "mathfem.h"
42 #include "masterdof.h"
43 #include "simpleslavedof.h"
44 #include "nonlocalbarrier.h"
45 #include "initialcondition.h"
46 #include "classfactory.h"
47 #include "outputmanager.h"
48 #include "crosssection.h"
49 #include "function.h"
50 #include "timer.h"
51 #include "remeshingcrit.h"
52 #include "dynamicinputrecord.h"
53 #include "engngm.h"
54 
55 #include <queue>
56 #include <set>
57 
58 #ifdef __OOFEG
59  #include "oofeggraphiccontext.h"
60 #endif
61 
62 #ifdef __PARALLEL_MODE
63  #include "parallel.h"
64  #include "problemcomm.h"
65  #include "communicator.h"
66  #include "datastream.h"
68 #endif
69 
70 namespace oofem {
71 
72 REGISTER_Mesher(Subdivision, MPT_SUBDIVISION);
73 
74 //#define __VERBOSE_PARALLEL
75 #define DEBUG_CHECK
76 //#define DEBUG_INFO
77 //#define DEBUG_SMOOTHING
78 
79 // QUICK_HACK enables comparison of sequential and parallel version when smoothing is applied;
80 // sequential nodes on parallel interface are marked as boundary (in input file);
81 // important: there is a limitation:
82 // single element is allowed to have only 2 (in 2D) or 3 (in 3D) nodes marked as boundary
83 // if quick hack is used, standard identification of boundary nodes during bisection is
84 // replaced by simplified version;
85 // also smoothing of nodes is driven by coordinates based ordering which should ensure
86 // that sequential nodes (on individual parallel partitions) are smoothed in the same
87 // order as parallel nodes
88 
89 //#define QUICK_HACK
90 
91 #ifdef __OOFEG
92  #define DRAW_IRREGULAR_NODES
93  #define DRAW_REMOTE_ELEMENTS
94  //#define DRAW_MESH_BEFORE_BISECTION
95  #define DRAW_MESH_AFTER_BISECTION
96  //#define DRAW_MESH_AFTER_EACH_BISECTION_LEVEL
97 #endif
98 
99 
100 
101 //#define NM // Nooru Mohamed
102 //#define BRAZIL_2D // splitting brazil test
103 //#define THREEPBT_3D // 3pbt
104 #define HEADEDSTUD // headed stud
105 
106 
107 
108 /* builds connectivity of top level node only however at any level of the subdivision */
109 
111 {
112  IntArray me(1), connElems;
113  int i, el;
114  me.at(1) = this->giveNumber();
115 
116  ct->giveNodeNeighbourList(connElems, me);
117  this->connectedElements.preallocate( connElems.giveSize() ); // estimate size of the list
118 
119  for ( i = 1; i <= connElems.giveSize(); i++ ) {
120  el = connElems.at(i);
121 
122 #ifdef __PARALLEL_MODE
123  if ( this->mesh->giveElement(el)->giveParallelMode() != Element_local ) {
124  continue;
125  }
126 #endif
127 
128  // use nonzero chunk, the estimated size may be not large enough
130  }
131 
132  return ( connectedElements.giveSize() );
133 }
134 
135 
137 {
138  int el, i;
140 
141  if ( this->isTerminal() ) {
142  node->insertConnectedElement( this->giveNumber() );
143  return;
144  }
145 
146  for ( el = 1; el <= this->giveChildren()->giveSize(); el++ ) {
147  elem = mesh->giveElement( this->giveChildren()->at(el) );
148  for ( i = 1; i <= elem->giveNodes()->giveSize(); i++ ) {
149  if ( elem->giveNode(i) == node->giveNumber() ) {
150  elem->buildTopLevelNodeConnectivity(node);
151  break;
152  }
153  }
154  }
155 }
156 
157 
158 #ifdef __PARALLEL_MODE
159 
160 int
162 {
163  IntArray me(1), connElems;
164  int i, el, cnt = 0;
165  me.at(1) = this->giveNumber();
166 
167  ct->giveNodeNeighbourList(connElems, me);
168  this->connectedElements.preallocate( connElems.giveSize() ); // estimate size of the list
169 
170  for ( i = 1; i <= connElems.giveSize(); i++ ) {
171  el = connElems.at(i);
172  if ( this->mesh->giveElement(el)->giveParallelMode() != Element_local ) {
173  continue;
174  }
175 
176  // use zero chunk, the array is large enough
177  this->connectedElements.insertSorted(el, 0);
178  cnt++;
179  }
180 
181  return ( cnt );
182 }
183 
184 
185 
186 void
188 {
189  int ie, num = this->giveNumber();
190  IntArray connNodes;
191  const IntArray *connElems = this->giveConnectedElements();
192 
193  for ( ie = 1; ie <= connElems->giveSize(); ie++ ) {
194  this->mesh->giveElement( connElems->at(ie) )->numberSharedEdges(num, connNodes);
195  }
196 }
197 
198 
199 void
201 {
202  int eIndex, jnode, jNode, eNum = this->mesh->giveNumberOfEdges();
204 
205  for ( jnode = 1; jnode <= 3; jnode++ ) {
206  jNode = nodes.at(jnode);
207  if ( jNode == iNode ) {
208  continue;
209  }
210 
211  if ( this->mesh->giveNode(jNode)->giveParallelMode() != DofManager_shared ) {
212  continue;
213  }
214 
215  // potential shared edge
216 
217  if ( connNodes.contains(jNode) ) {
218  continue; // edge already processed (from iNode)
219  }
220 
221  connNodes.followedBy(jNode, 10);
222 
223  eIndex = giveEdgeIndex(iNode, jNode);
224 
225  if ( shared_edges.giveSize() ) {
226  if ( shared_edges.at(eIndex) ) {
227  continue; // edge already processed (from jNode)
228  }
229  }
230 
231  // create edge (even if it might not be shared)
232  _edge = new Subdivision :: RS_SharedEdge(this->mesh);
233  _edge->setEdgeNodes(iNode, jNode);
234 
235  this->mesh->addEdge(_edge);
236  eNum++;
237 
238  // get the tentative partitions
239  // they must be confirmed via mutual communication
241  if ( _edge->giveSharedPartitions(partitions) ) {
242  _edge->setPartitions(partitions);
243 
244  // I do rely on the fact that in 2D shared edge can be incident only to one element !!!
245  if ( !shared_edges.giveSize() ) {
246  makeSharedEdges();
247  }
248 
249  shared_edges.at(eIndex) = eNum;
250  }
251  }
252 }
253 
254 
255 void
257 {
258  int ie, elems, eIndex, jnode, jNode, eNum = this->mesh->giveNumberOfEdges();
261 
262  for ( jnode = 1; jnode <= 4; jnode++ ) {
263  jNode = nodes.at(jnode);
264  if ( jNode == iNode ) {
265  continue;
266  }
267 
268  if ( this->mesh->giveNode(jNode)->giveParallelMode() != DofManager_shared ) {
269  continue;
270  }
271 
272  // potential shared edge
273 
274  if ( connNodes.contains(jNode) ) {
275  continue; // edge already processed (from iNode)
276  }
277 
278  connNodes.followedBy(jNode, 20);
279 
280  eIndex = giveEdgeIndex(iNode, jNode);
281 
282  if ( shared_edges.giveSize() ) {
283  if ( shared_edges.at(eIndex) ) {
284  continue; // edge already processed (from jNode)
285  }
286  }
287 
288  // create edge (even if it might not be shared)
289  _edge = new Subdivision :: RS_SharedEdge(this->mesh);
290  _edge->setEdgeNodes(iNode, jNode);
291 
292  this->mesh->addEdge(_edge);
293  ++eNum;
294 
295  // get the tentative partitions
296  // they must be confirmed via mutual communication
298  if ( _edge->giveSharedPartitions(partitions) ) {
299  _edge->setPartitions(partitions);
300  if ( !this->giveSharedEdges()->giveSize() ) {
301  this->makeSharedEdges();
302  }
303 
304  shared_edges.at(eIndex) = eNum;
305 
306  // put edge on all elements sharing it
307  const IntArray *iElems, *jElems;
308 
309  iElems = mesh->giveNode(iNode)->giveConnectedElements();
310  jElems = mesh->giveNode(jNode)->giveConnectedElements();
311 
312  IntArray common;
313  if ( iElems->giveSize() <= jElems->giveSize() ) {
314  common.preallocate( iElems->giveSize() );
315  } else {
316  common.preallocate( jElems->giveSize() );
317  }
318 
319  // I do rely on the fact that the arrays are ordered !!!
320  // I am using zero chunk because common is large enough
321  elems = iElems->findCommonValuesSorted(* jElems, common, 0);
322  #ifdef DEBUG_CHECK
323  if ( !elems ) {
324  OOFEM_ERROR("potentionally shared edge %d is not shared by common elements", eNum);
325  }
326 
327  #endif
328  for ( ie = 1; ie <= elems; ie++ ) {
329  elem = mesh->giveElement( common.at(ie) );
330  if ( elem == this ) {
331  continue;
332  }
333 
334  eIndex = elem->giveEdgeIndex(iNode, jNode);
335  if ( !elem->giveSharedEdges()->giveSize() ) {
336  elem->makeSharedEdges();
337  }
338 
339  elem->setSharedEdge(eIndex, eNum);
340  }
341  }
342  }
343 }
344 
345 
346 // CAUTION: gives the tentative partitions;
347 // they must be confirmed via mutual communication
348 
349 int
351 {
352  int count = 0;
353  int myrank = this->mesh->giveSubdivision()->giveRank();
354  const IntArray *iPartitions = mesh->giveNode(iNode)->givePartitions();
355  const IntArray *jPartitions = mesh->giveNode(jNode)->givePartitions();
356 
357  // find common partitions;
358  // this could be done more efficiently (using method findCommonValuesSorted)
359  // if array iPartitions and jPartitions were sorted;
360  // however they are not sorted in current version
361  int ipsize = iPartitions->giveSize();
362  for ( int i = 1; i <= ipsize; i++ ) {
363  if ( jPartitions->contains( iPartitions->at(i) ) && ( iPartitions->at(i) != myrank ) ) {
364  partitions.followedBy(iPartitions->at(i), 2);
365  count++;
366  }
367  }
368 
369  return count;
370 }
371 
372 #endif
373 
374 
375 double
377 {
378  double answer = 0.0;
379  for ( int n: nodes ) {
380  answer += mesh->giveNode( n )->giveRequiredDensity();
381  }
382 
383  answer = answer / nodes.giveSize();
384  return answer;
385 }
386 
387 
388 int
390 {
391  RS_Element *elem = this;
392  while ( elem->giveNumber() != elem->giveParent() ) {
393  elem = mesh->giveElement( elem->giveParent() );
394  }
395 
396  return ( elem->giveNumber() );
397 }
398 
399 
401  Subdivision :: RS_Element(number, mesh, parent, nodes)
402 {
407 
408 #ifdef DEBUG_CHECK
409  if ( nodes.findFirstIndexOf(0) ) {
410  OOFEM_ERROR("0 node of element %d", this->number);
411  }
412 
413 #endif
414 }
415 
416 
418  Subdivision :: RS_Element(number, mesh, parent, nodes)
419 {
422  side_leIndex.resize(4);
423  side_leIndex.zero();
426 
427 #ifdef DEBUG_CHECK
428  if ( nodes.findFirstIndexOf(0) ) {
429  OOFEM_ERROR("0 node of element %d", this->number);
430  }
431 
432 #endif
433 }
434 
435 
436 int
438 {
439  if ( !this->leIndex ) { // prevent multiple evaluation
440  double elength1, elength2, elength3;
441 
442  elength1 = mesh->giveNode( this->nodes.at(1) )->giveCoordinates()->distance_square( * ( mesh->giveNode( this->nodes.at(2) )->giveCoordinates() ) );
443  elength2 = mesh->giveNode( this->nodes.at(2) )->giveCoordinates()->distance_square( * ( mesh->giveNode( this->nodes.at(3) )->giveCoordinates() ) );
444  elength3 = mesh->giveNode( this->nodes.at(3) )->giveCoordinates()->distance_square( * ( mesh->giveNode( this->nodes.at(1) )->giveCoordinates() ) );
445 
446  leIndex = 1;
447  if ( elength2 > elength1 ) {
448  leIndex = 2;
449  if ( elength3 > elength2 ) {
450  leIndex = 3;
451  }
452  } else if ( elength3 > elength1 ) {
453  leIndex = 3;
454  }
455  }
456 
457  return ( leIndex );
458 }
459 
460 
461 int
463 {
464  if ( !this->leIndex ) { // prevent multiple evaluation
465  // IMPORTANT: ambiguity must be handled !!!
466 
467  // in order to achieve that shared tetra faces are subdivided in compatible way,
468  // always the same longest edge must be identified;
469  // this may be problem if there are two or more edges of the "same" length
470  // then the longest is dependent on the order in which edges are processed;
471  // also in order to get always the same distance of two nodes, the nodes should be processes in the same order;
472 
473  // however it is not clear whether one may rely that the same distance is always obtained for the same pair of nodes
474  // (in proper order) especially if processed on different processors
475  // in parallel environment some additional communication may be needed !!! HUHU
476 
477  // if the second largest edge is not equal to the largest edge and if the largest and second largest edges
478  // are opposite edges, no action is theoretically required;
479  // but because it is not clear whether different result would be obtained if edge lengths are evaluated
480  // using node ordering, this possibility (of no action) is not considered
481 
482  double elength [ 6 ], maxlength, maxlen [ 4 ];
483  int i, j, side, l_nd [ 4 ], ind [ 6 ];
484  // array ed_side contains face numbers shared by the edge (indexing from 0)
485  int nd [ 4 ] = {
486  0, 1, 2, 3
487  }, ed_side [ 6 ] [ 2 ] = { { 0, 1 }, { 0, 2 }, { 0, 3 }, { 3, 1 }, { 1, 2 }, { 2, 3 } };
488  // array ed contains edge numbers (1 to 6) connecting the node with each node (1 to 4)
489  // (indexing from 1, zero means, that node is not connected by an edge with itself)
490  int ed [ 4 ] [ 4 ] = { { 0, 1, 3, 4 }, { 1, 0, 2, 5 }, { 3, 2, 0, 6 }, { 4, 5, 6, 0 } };
491 #ifdef __PARALLEL_MODE
492  int g_nd [ 4 ];
493 #endif
494  bool swap;
495 
496  // sort node ids
497  for ( i = 0; i < 4; i++ ) {
498  l_nd [ i ] = nodes.at(i + 1);
499 #ifdef __PARALLEL_MODE
500  g_nd [ i ] = mesh->giveNode(l_nd [ i ])->giveGlobalNumber();
501 #endif
502  }
503 
504  swap = true;
505  while ( swap ) {
506  swap = false;
507  for ( i = 0; i < 3; i++ ) {
508 #ifdef __PARALLEL_MODE
509  if ( g_nd [ i ] > g_nd [ i + 1 ] ) {
510 #else
511  if ( l_nd [ i ] > l_nd [ i + 1 ] ) {
512 #endif
513  int tmp;
514 #ifdef __PARALLEL_MODE
515  tmp = g_nd [ i ];
516  g_nd [ i ] = g_nd [ i + 1 ];
517  g_nd [ i + 1 ] = tmp;
518 #endif
519  tmp = l_nd [ i ];
520  l_nd [ i ] = l_nd [ i + 1 ];
521  l_nd [ i + 1 ] = tmp;
522  tmp = nd [ i ];
523  nd [ i ] = nd [ i + 1 ];
524  nd [ i + 1 ] = tmp;
525  swap = true;
526  }
527  }
528  }
529 
530  // assign edge indices from smallest node id !!!
531  ind [ 0 ] = ed [ nd [ 0 ] ] [ nd [ 1 ] ];
532  ind [ 1 ] = ed [ nd [ 0 ] ] [ nd [ 2 ] ];
533  ind [ 2 ] = ed [ nd [ 0 ] ] [ nd [ 3 ] ];
534  ind [ 3 ] = ed [ nd [ 1 ] ] [ nd [ 2 ] ];
535  ind [ 4 ] = ed [ nd [ 1 ] ] [ nd [ 3 ] ];
536  ind [ 5 ] = ed [ nd [ 2 ] ] [ nd [ 3 ] ];
537 
538  // evaluate edge lengths from smallest node id !!!
539  elength [ 0 ] = mesh->giveNode(l_nd [ 0 ])->giveCoordinates()->distance_square( * ( mesh->giveNode(l_nd [ 1 ])->giveCoordinates() ) );
540  elength [ 1 ] = mesh->giveNode(l_nd [ 0 ])->giveCoordinates()->distance_square( * ( mesh->giveNode(l_nd [ 2 ])->giveCoordinates() ) );
541  elength [ 2 ] = mesh->giveNode(l_nd [ 0 ])->giveCoordinates()->distance_square( * ( mesh->giveNode(l_nd [ 3 ])->giveCoordinates() ) );
542  elength [ 3 ] = mesh->giveNode(l_nd [ 1 ])->giveCoordinates()->distance_square( * ( mesh->giveNode(l_nd [ 2 ])->giveCoordinates() ) );
543  elength [ 4 ] = mesh->giveNode(l_nd [ 1 ])->giveCoordinates()->distance_square( * ( mesh->giveNode(l_nd [ 3 ])->giveCoordinates() ) );
544  elength [ 5 ] = mesh->giveNode(l_nd [ 2 ])->giveCoordinates()->distance_square( * ( mesh->giveNode(l_nd [ 3 ])->giveCoordinates() ) );
545 
546  // get absolutely largest edge and the largest edge on individual sides
547  // process edges in the same order in which their length was evaluated !!!
548  maxlength = maxlen [ 0 ] = maxlen [ 1 ] = maxlen [ 2 ] = maxlen [ 3 ] = 0.0;
549  for ( i = 0; i < 6; i++ ) {
550  if ( elength [ i ] > maxlength ) {
551  maxlength = elength [ i ];
552  leIndex = ind [ i ];
553  for ( j = 0; j < 2; j++ ) {
554  side = ed_side [ ind [ i ] - 1 ] [ j ];
555  maxlen [ side ] = elength [ i ];
556  side_leIndex.at(side + 1) = ind [ i ];
557  }
558  } else {
559  for ( j = 0; j < 2; j++ ) {
560  side = ed_side [ ind [ i ] - 1 ] [ j ];
561  if ( elength [ i ] > maxlen [ side ] ) {
562  maxlen [ side ] = elength [ i ];
563  side_leIndex.at(side + 1) = ind [ i ];
564  }
565  }
566  }
567  }
568  }
569 
570  return ( leIndex );
571 }
572 
573 
574 int
576 {
577  /* returns zero, if triangle does not have iNode or jNode */
578  int in = 0, jn = 0;
579 
580  int k;
581  for ( k = 1; k <= 3; k++ ) {
582  if ( nodes.at(k) == iNode ) {
583  in = k;
584  }
585 
586  if ( nodes.at(k) == jNode ) {
587  jn = k;
588  }
589  }
590 
591  if ( in == 0 || jn == 0 ) {
592  OOFEM_ERROR( "there is no edge connecting %d and %d on element %d",
593  iNode, jNode, this->giveNumber() );
594  return 0;
595  }
596 
597  if ( in < jn ) {
598  if ( jn == in + 1 ) {
599  return ( in );
600  }
601 
602  return ( 3 );
603  } else {
604  if ( in == jn + 1 ) {
605  return ( jn );
606  }
607 
608  return ( 3 );
609  }
610 }
611 
612 
613 int
615 {
616  /* returns zero, if tetra does not have iNode or jNode */
617  int in = 0, jn = 0;
618 
619  int k;
620  for ( k = 1; k <= 4; k++ ) {
621  if ( nodes.at(k) == iNode ) {
622  in = k;
623  }
624 
625  if ( nodes.at(k) == jNode ) {
626  jn = k;
627  }
628  }
629 
630  if ( in == 0 || jn == 0 ) {
631  OOFEM_ERROR( "there is no edge connecting %d and %d on element %d",
632  iNode, jNode, this->giveNumber() );
633  return 0;
634  }
635 
636  if ( in < jn ) {
637  if ( jn == 4 ) {
638  return ( in + 3 );
639  }
640 
641  if ( jn == in + 1 ) {
642  return ( in );
643  }
644 
645  return ( 3 );
646  } else {
647  if ( in == 4 ) {
648  return ( jn + 3 );
649  }
650 
651  if ( in == jn + 1 ) {
652  return ( jn );
653  }
654 
655  return ( 3 );
656  }
657 }
658 
659 
660 void
662 {
663  /* this is symbolic bisection - no new elements are added, only irregular nodes are introduced */
664  int inode, jnode, iNode, jNode, iNum, eInd;
665  double density;
666  bool boundary = false;
667  FloatArray coords;
669 #ifdef __PARALLEL_MODE
671 #endif
672 
673  if ( !irregular_nodes.at(leIndex) ) {
674  RS_IrregularNode *irregular;
675  // irregular on the longest edge does not exist
676  // introduce new irregular node on the longest edge
677  inode = leIndex;
678  jnode = ( inode < 3 ) ? inode + 1 : 1;
679 
680  iNode = nodes.at(inode);
681  jNode = nodes.at(jnode);
682  // compute coordinates of new irregular
683  coords = * ( mesh->giveNode(iNode)->giveCoordinates() );
684  coords.add( * mesh->giveNode(jNode)->giveCoordinates() );
685  coords.times(0.5);
686  // compute required density of a new node
687  density = 0.5 * ( mesh->giveNode(iNode)->giveRequiredDensity() +
688  mesh->giveNode(jNode)->giveRequiredDensity() );
689  // create new irregular
690  iNum = mesh->giveNumberOfNodes() + 1;
691  irregular = new Subdivision :: RS_IrregularNode(iNum, mesh, 0, coords, density, boundary);
692  mesh->addNode(irregular);
693  // add irregular to receiver
694  this->irregular_nodes.at(leIndex) = iNum;
695 
696 #ifdef __OOFEG
697  #ifdef DRAW_IRREGULAR_NODES
698  irregular->drawGeometry();
699  #endif
700 #endif
701 
702 #ifdef __PARALLEL_MODE
703 #ifdef __VERBOSE_PARALLEL
704  OOFEM_LOG_INFO("[%d] RS_Triangle::bisecting %d nodes %d %d %d, leIndex %d, new irregular %d\n", mesh->giveSubdivision()->giveRank(), this->number, nodes.at(1), nodes.at(2), nodes.at(3), leIndex, iNum);
705 #endif
706 #endif
707 
708 #ifdef QUICK_HACK
709  if ( mesh->giveNode( nodes.at(1) )->isBoundary() && mesh->giveNode( nodes.at(2) )->isBoundary() && mesh->giveNode( nodes.at(3) )->isBoundary() ) {
710  OOFEM_ERROR( "quick hack not applicable due to element %d", this->giveNumber() );
711  }
712 
713  if ( mesh->giveNode(iNode)->isBoundary() && mesh->giveNode(jNode)->isBoundary() ) {
714  boundary = true;
715  }
716 
717 #else
718  // check whether new node is boundary
719  if ( this->neghbours_base_elements.at(leIndex) ) {
720  Domain *dorig = mesh->giveSubdivision()->giveDomain();
721  // I rely on tha fact that nodes on intermaterial interface are marked as boundary
722  // however this might not be true
723  // therefore in smoothing the boundary flag is setuped again if not set boundary from here
724  if ( mesh->giveNode(iNode)->isBoundary() || mesh->giveNode(jNode)->isBoundary() ) {
725  if ( dorig->giveElement( this->giveTopParent() )->giveRegionNumber() != dorig->giveElement( mesh->giveElement( this->neghbours_base_elements.at(leIndex) )->giveTopParent() )->giveRegionNumber() ) {
726  boundary = true;
727  }
728  }
729  } else {
730  boundary = true;
731  }
732 
733 #endif
734 
735  if ( boundary ) {
736  irregular->setBoundary(true);
737  }
738 
739  if ( this->neghbours_base_elements.at(leIndex) ) {
740  // add irregular to neighbour
742  eInd = elem->giveEdgeIndex(iNode, jNode);
743  elem->setIrregular(eInd, iNum);
744 
745  if ( !elem->giveQueueFlag() ) {
746  // add neighbour to list of elements for subdivision
747  subdivqueue.push( this->neghbours_base_elements.at(leIndex) );
748  elem->setQueueFlag(true);
749  }
750  }
751 
752 #ifdef __PARALLEL_MODE
753  else {
754  // check if there are (potentionally) shared edges
755  if ( shared_edges.giveSize() ) {
756  // check if the edge is (really) shared
757  if ( shared_edges.at(leIndex) ) {
758  edge = mesh->giveEdge( shared_edges.at(leIndex) );
759 
760  #ifdef DEBUG_CHECK
761  if ( !edge->givePartitions()->giveSize() ) {
762  OOFEM_ERROR( "unshared edge %d of element %d is marked as shared",
763  shared_edges.at(leIndex), this->giveNumber() );
764  }
765 
766  #endif
767 
768  // new node is on shared interpartition boundary
770  // partitions are inherited from shared edge
771  irregular->setPartitions( * ( edge->givePartitions() ) );
772  irregular->setEdgeNodes(iNode, jNode);
773  // put its number into queue of shared irregulars that is later used to inform remote partitions about this fact
774  sharedIrregularsQueue.push_back(iNum);
775  #ifdef __VERBOSE_PARALLEL
776  OOFEM_LOG_INFO("RS_Triangle::bisect: Shared irregular detected, number %d nodes %d %d [%d %d], elem %d\n", iNum, iNode, jNode, mesh->giveNode(iNode)->giveGlobalNumber(), mesh->giveNode(jNode)->giveGlobalNumber(), this->number);
777  #endif
778  }
779  }
780  }
781 #endif
782  }
783 
784  this->setQueueFlag(false);
785 }
786 
787 
788 void
789 Subdivision :: RS_Tetra :: bisect(std :: queue< int > &subdivqueue, std :: list< int > &sharedIrregularsQueue)
790 {
791  /* this is symbolic bisection - no new elements are added, only irregular nodes are introduced */
792  int i, j, inode, jnode, iNode, jNode, ngb, eIndex, eInd, reg, iNum, side, cnt = 0, elems;
793  // array ed_side contains face numbers NOT shared by the edge (indexing from 1)
794  int ed_side [ 6 ] [ 2 ] = { { 3, 4 }, { 4, 2 }, { 2, 3 }, { 1, 3 }, { 1, 4 }, { 1, 2 } }, ed [ 4 ] = {
795  0, 0, 0, 0
796  }, opp_ed [ 6 ] = {
797  6, 4, 5, 2, 3, 1
798  };
799  // array side_ed contains edge numbers bounding faces (indexing from 1)
800  int side_ed [ 4 ] [ 3 ] = { { 1, 2, 3 }, { 1, 5, 4 }, { 2, 6, 5 }, { 3, 4, 6 } };
801  double density;
802  bool shared, boundary, iboundary, jboundary, opposite = false;
803  Subdivision :: RS_Element *elem1, *elem2, *elem;
804  FloatArray coords;
805  Domain *dorig;
806  RS_IrregularNode *irregular;
807  const IntArray *iElems, *jElems;
808 #ifdef __PARALLEL_MODE
810 #endif
811 
812  dorig = mesh->giveSubdivision()->giveDomain();
813  reg = dorig->giveElement( this->giveTopParent() )->giveRegionNumber();
814 
815  // first resolve whether there will be inserted irregular on the edge opposite to longest edge
816  // this will happen if the opposite edge is longest for a side not shared by the longest edge
817  // and simultaneously there is already an irregular on that side
818  if ( !irregular_nodes.at(opp_ed [ leIndex - 1 ]) ) {
819  for ( i = 0; i < 2; i++ ) {
820  // get side not incident to leIndex edge
821  side = ed_side [ leIndex - 1 ] [ i ];
822  // check whether the longest edge of that side is opposite edge
823  if ( side_leIndex.at(side) != opp_ed [ leIndex - 1 ] ) {
824  continue;
825  }
826 
827  for ( j = 0; j < 3; j++ ) {
828  // check for irregulars on that side
829  if ( irregular_nodes.at(side_ed [ side - 1 ] [ j ]) ) {
830  opposite = true;
831  irregular_nodes.at(opp_ed [ leIndex - 1 ]) = -1; // fictitious irregular
832  break;
833  }
834  }
835 
836  if ( opposite ) {
837  break;
838  }
839  }
840  }
841 
842  // insert irregular on absolutely longest edge first;
843  // this controls the primary bisection;
844  // if there is an irregular (including a fictitious one) on side not shared by longest edge
845  // insert irregular on longest edge of that side;
846  ed [ cnt++ ] = leIndex;
847  for ( i = 0; i < 2; i++ ) {
848  // get side not incident to leIndex edge
849  side = ed_side [ leIndex - 1 ] [ i ];
850  for ( j = 0; j < 3; j++ ) {
851  // check for irregulars on that side
852  if ( irregular_nodes.at(side_ed [ side - 1 ] [ j ]) ) {
853  ed [ cnt++ ] = side_leIndex.at(side);
854  break;
855  }
856  }
857  }
858 
859  // remove fictitious irregular
860  if ( opposite ) {
861  irregular_nodes.at(opp_ed [ leIndex - 1 ]) = 0;
862  }
863 
864  // check for the case when longest edge of the sides not shared by the absolutely longest edge is the same
865  // (it is opposite to absolutely longest edge)
866  if ( cnt == 3 ) {
867  if ( ed [ 1 ] == ed [ 2 ] ) {
868  cnt = 2;
869 #ifdef DEBUG_CHECK
870  if ( ed [ 1 ] != opp_ed [ leIndex - 1 ] ) {
871  OOFEM_ERROR( "unexpected situation on element %d", this->giveNumber() );
872  }
873 
874 #endif
875  }
876  }
877 
878 #ifdef QUICK_HACK
879  if ( mesh->giveNode( nodes.at(1) )->isBoundary() && mesh->giveNode( nodes.at(2) )->isBoundary() &&
880  mesh->giveNode( nodes.at(3) )->isBoundary() && mesh->giveNode( nodes.at(4) )->isBoundary() ) {
881  OOFEM_ERROR( "quick hack not applicable due to element %d", this->giveNumber() );
882  }
883 
884 #endif
885 
886  // insert all relevant irregulars
887  for ( i = 0; i < cnt; i++ ) {
888  eIndex = ed [ i ];
889  if ( !irregular_nodes.at(eIndex) ) {
890  // irregular on this edge does not exist;
891  // introduce new irregular node on this edge on this element and on all local elements sharing that edge;
892  // if the edge is local, neigbours are processed;
893  // if the edge is shared, elements sharing simultaneously both end nodes are processed;
894  if ( eIndex <= 3 ) {
895  inode = eIndex;
896  jnode = ( eIndex < 3 ) ? inode + 1 : 1;
897  ngb = 1;
898  } else {
899  inode = eIndex - 3;
900  jnode = 4;
901  ngb = inode + 1;
902  }
903 
904  iNode = nodes.at(inode);
905  jNode = nodes.at(jnode);
906  // compute coordinates of new irregular
907  coords = * ( mesh->giveNode(iNode)->giveCoordinates() );
908  coords.add( * mesh->giveNode(jNode)->giveCoordinates() );
909  coords.times(0.5);
910 #ifdef HEADEDSTUD
911  double dist, rad, rate;
912  FloatArray *c;
913 
914  c = mesh->giveNode(iNode)->giveCoordinates();
915  dist = c->at(1) * c->at(1) + c->at(3) * c->at(3);
916  if ( c->at(2) > 69.9999999 ) {
917  rad = 7.0;
918  } else if ( c->at(2) < 64.5000001 ) {
919  rad = 18.0;
920  } else {
921  rad = 18.0 - 11.0 / 5.5 * ( c->at(2) - 64.5 );
922  }
923 
924  if ( fabs(dist - rad * rad) < 0.01 ) { // be very tolerant (geometry is not precise)
925  c = mesh->giveNode(jNode)->giveCoordinates();
926  dist = c->at(1) * c->at(1) + c->at(3) * c->at(3);
927  if ( c->at(2) > 69.9999999 ) {
928  rad = 7.0;
929  } else if ( c->at(2) < 64.5000001 ) {
930  rad = 18.0;
931  } else {
932  rad = 18.0 - 11.0 / 5.5 * ( c->at(2) - 64.5 );
933  }
934 
935  if ( fabs(dist - rad * rad) < 0.01 ) { // be very tolerant (geometry is not precise)
936  dist = coords.at(1) * coords.at(1) + coords.at(3) * coords.at(3);
937  if ( coords.at(2) > 69.9999999 ) {
938  rad = 7.0;
939  } else if ( coords.at(2) < 64.5000001 ) {
940  rad = 18.0;
941  } else {
942  rad = 18.0 - 11.0 / 5.5 * ( coords.at(2) - 64.5 );
943  }
944 
945  rate = rad / sqrt(dist);
946  coords.at(1) *= rate;
947  coords.at(3) *= rate;
948  }
949  }
950 
951 #endif
952  // compute required density of a new node
953  density = 0.5 * ( mesh->giveNode(iNode)->giveRequiredDensity() +
954  mesh->giveNode(jNode)->giveRequiredDensity() );
955  // create new irregular
956  iNum = mesh->giveNumberOfNodes() + 1;
957  irregular = new Subdivision :: RS_IrregularNode(iNum, mesh, 0, coords, density, false);
958  mesh->addNode(irregular);
959  // add irregular to receiver
960  this->irregular_nodes.at(eIndex) = iNum;
961 
962 #ifdef __OOFEG
963  #ifdef DRAW_IRREGULAR_NODES
964  irregular->drawGeometry();
965  #endif
966 #endif
967 
968 #ifdef DEBUG_INFO
969  #ifdef __PARALLEL_MODE
970  // do not print global numbers of elements because they are not available (they are assigned at once after bisection);
971  // do not print global numbers of irregulars as these may not be available yet
972  OOFEM_LOG_INFO( "[%d] Irregular %d added on %d [%d] (edge %d, nodes %d %d [%d %d], nds %d %d %d %d [%d %d %d %d], ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
973  mesh->giveSubdivision()->giveRank(), iNum, this->number, this->giveGlobalNumber(), eIndex, iNode, jNode,
975  nodes.at(1), nodes.at(2), nodes.at(3), nodes.at(4),
982  #else
983  OOFEM_LOG_INFO( "Irregular %d added on %d (edge %d, nodes %d %d, nds %d %d %d %d, ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
984  iNum, this->number, eIndex, iNode, jNode,
985  nodes.at(1), nodes.at(2), nodes.at(3), nodes.at(4),
990  #endif
991 #endif
992 
993  shared = boundary = false;
994 
995 #ifdef __PARALLEL_MODE
996  // check if there are (potentionally) shared edges
997  if ( shared_edges.giveSize() ) {
998  // check if the edge is (really) shared
999  if ( shared_edges.at(eIndex) ) {
1000  edge = mesh->giveEdge( shared_edges.at(eIndex) );
1001 
1002  #ifdef DEBUG_CHECK
1003  if ( !edge->givePartitions()->giveSize() ) {
1004  OOFEM_ERROR( "unshared edge %d of element %d is marked as shared",
1005  shared_edges.at(eIndex), this->giveNumber() );
1006  }
1007 
1008  #endif
1009 
1010  shared = boundary = true;
1011  // new node is on shared interpartition boundary
1012  irregular->setParallelMode(DofManager_shared);
1013  irregular->setPartitions( * ( edge->givePartitions() ) );
1014  irregular->setEdgeNodes(iNode, jNode);
1015  // put its number into queue of shared irregulars that is later used to inform remote partitions about this fact
1016  sharedIrregularsQueue.push_back(iNum);
1017  #ifdef __VERBOSE_PARALLEL
1018  OOFEM_LOG_INFO("RS_Tetra::bisect: Shared irregular detected, number %d nodes %d %d [%d %d], elem %d\n", iNum, iNode, jNode, mesh->giveNode(iNode)->giveGlobalNumber(), mesh->giveNode(jNode)->giveGlobalNumber(), this->number);
1019  #endif
1020 
1021  iElems = mesh->giveNode(iNode)->giveConnectedElements();
1022  jElems = mesh->giveNode(jNode)->giveConnectedElements();
1023 
1024  IntArray common;
1025  if ( iElems->giveSize() <= jElems->giveSize() ) {
1026  common.preallocate( iElems->giveSize() );
1027  } else {
1028  common.preallocate( jElems->giveSize() );
1029  }
1030 
1031  // I do rely on the fact that the arrays are ordered !!!
1032  // I am using zero chunk because common is large enough
1033  elems = iElems->findCommonValuesSorted(* jElems, common, 0);
1034  #ifdef DEBUG_CHECK
1035  if ( !elems ) {
1036  OOFEM_ERROR( "shared edge %d is not shared by common elements",
1037  shared_edges.at(eIndex) );
1038  }
1039 
1040  #endif
1041  // after subdivision there will be twice as much of connected local elements
1042  irregular->preallocateConnectedElements(elems * 2);
1043  // put the new node on appropriate edge of all local elements (except "this") sharing both nodes
1044  for ( j = 1; j <= elems; j++ ) {
1045  elem = mesh->giveElement( common.at(j) );
1046  if ( elem == this ) {
1047  continue;
1048  }
1049 
1050  #ifdef DEBUG_CHECK
1051  if ( !elem->giveSharedEdges()->giveSize() ) {
1052  OOFEM_ERROR( "element %d incident to shared edge %d does not have shared edges",
1053  common.at(j), shared_edges.at(eIndex) );
1054  }
1055 
1056  #endif
1057 
1058  eInd = elem->giveEdgeIndex(iNode, jNode);
1059  elem->setIrregular(eInd, iNum);
1060 
1061  if ( !elem->giveQueueFlag() ) {
1062  // add elem to list of elements for subdivision
1063  subdivqueue.push( elem->giveNumber() );
1064  elem->setQueueFlag(true);
1065  }
1066 
1067  #ifdef DEBUG_INFO
1068  #ifdef __PARALLEL_MODE
1069  // do not print global numbers of elements because they are not available (they are assigned at once after bisection);
1070  // do not print global numbers of irregulars as these may not be available yet
1071  OOFEM_LOG_INFO( "[%d] Irregular %d added on %d [%d] (edge %d, nodes %d %d [%d %d], nds %d %d %d %d [%d %d %d %d], ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
1072  mesh->giveSubdivision()->giveRank(), iNum,
1073  elem->giveNumber(), this->giveGlobalNumber(), eInd, iNode, jNode,
1075  elem->giveNode(1), elem->giveNode(2), elem->giveNode(3), elem->giveNode(4),
1076  mesh->giveNode( elem->giveNode(1) )->giveGlobalNumber(),
1077  mesh->giveNode( elem->giveNode(2) )->giveGlobalNumber(),
1078  mesh->giveNode( elem->giveNode(3) )->giveGlobalNumber(),
1079  mesh->giveNode( elem->giveNode(4) )->giveGlobalNumber(),
1080  elem->giveNeighbor(1), elem->giveNeighbor(2), elem->giveNeighbor(3), elem->giveNeighbor(4),
1081  elem->giveIrregular(1), elem->giveIrregular(2), elem->giveIrregular(3),
1082  elem->giveIrregular(4), elem->giveIrregular(5), elem->giveIrregular(6) );
1083  #else
1084  OOFEM_LOG_INFO( "Irregular %d added on %d (edge %d, nodes %d %d, nds %d %d %d %d, ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
1085  iNum, elem->giveNumber(), eInd, iNode, jNode,
1086  elem->giveNode(1), elem->giveNode(2), elem->giveNode(3), elem->giveNode(4),
1087  elem->giveNeighbor(1), elem->giveNeighbor(2), elem->giveNeighbor(3), elem->giveNeighbor(4),
1088  elem->giveIrregular(1), elem->giveIrregular(2), elem->giveIrregular(3),
1089  elem->giveIrregular(4), elem->giveIrregular(5), elem->giveIrregular(6) );
1090  #endif
1091  #endif
1092  }
1093  }
1094  }
1095 
1096 #endif
1097 
1098  if ( !shared ) {
1099  iboundary = mesh->giveNode(iNode)->isBoundary();
1100  jboundary = mesh->giveNode(jNode)->isBoundary();
1101 #ifdef QUICK_HACK
1102  if ( iboundary == true && jboundary == true ) {
1103  boundary = true;
1104  }
1105 
1106 #endif
1107 
1108  // traverse neighbours
1109  elem1 = this;
1110  elem2 = NULL;
1111  while ( elem1->giveNeighbor(ngb) ) {
1112  elem2 = mesh->giveElement( elem1->giveNeighbor(ngb) );
1113  if ( elem2 == this ) {
1114  break;
1115  }
1116 
1117  eInd = elem2->giveEdgeIndex(iNode, jNode);
1118  elem2->setIrregular(eInd, iNum);
1119 
1120  if ( !elem2->giveQueueFlag() ) {
1121  // add neighbour to list of elements for subdivision
1122  subdivqueue.push( elem2->giveNumber() );
1123  elem2->setQueueFlag(true);
1124  }
1125 
1126 #ifndef QUICK_HACK
1127  if ( !boundary ) {
1128  // I rely on the fact that nodes on intermaterial interface are marked as boundary
1129  // however this might not be true
1130  // therefore in smoothing the boundary flag is setuped again if not set boundary from here
1131  if ( iboundary == true || jboundary == true ) {
1132  if ( dorig->giveElement( elem2->giveTopParent() )->giveRegionNumber() != reg ) {
1133  boundary = true;
1134  }
1135  }
1136  }
1137 
1138 #endif
1139 
1140 #ifdef DEBUG_INFO
1141  #ifdef __PARALLEL_MODE
1142  // do not print global numbers of elements because they are not available (they are assigned at once after bisection);
1143  // do not print global numbers of irregulars as these may not be available yet
1144  OOFEM_LOG_INFO( "[%d] Irregular %d added on %d [%d] (edge %d, nodes %d %d [%d %d], nds %d %d %d %d [%d %d %d %d], ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
1145  mesh->giveSubdivision()->giveRank(), iNum,
1146  elem2->giveNumber(), this->giveGlobalNumber(), eInd, iNode, jNode,
1148  elem2->giveNode(1), elem2->giveNode(2), elem2->giveNode(3), elem2->giveNode(4),
1149  mesh->giveNode( elem2->giveNode(1) )->giveGlobalNumber(),
1150  mesh->giveNode( elem2->giveNode(2) )->giveGlobalNumber(),
1151  mesh->giveNode( elem2->giveNode(3) )->giveGlobalNumber(),
1152  mesh->giveNode( elem2->giveNode(4) )->giveGlobalNumber(),
1153  elem2->giveNeighbor(1), elem2->giveNeighbor(2), elem2->giveNeighbor(3), elem2->giveNeighbor(4),
1154  elem2->giveIrregular(1), elem2->giveIrregular(2), elem2->giveIrregular(3),
1155  elem2->giveIrregular(4), elem2->giveIrregular(5), elem2->giveIrregular(6) );
1156  #else
1157  OOFEM_LOG_INFO( "Irregular %d added on %d (edge %d, nodes %d %d, nds %d %d %d %d, ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
1158  iNum, elem2->giveNumber(), eInd, iNode, jNode,
1159  elem2->giveNode(1), elem2->giveNode(2), elem2->giveNode(3), elem2->giveNode(4),
1160  elem2->giveNeighbor(1), elem2->giveNeighbor(2), elem2->giveNeighbor(3), elem2->giveNeighbor(4),
1161  elem2->giveIrregular(1), elem2->giveIrregular(2), elem2->giveIrregular(3),
1162  elem2->giveIrregular(4), elem2->giveIrregular(5), elem2->giveIrregular(6) );
1163  #endif
1164 #endif
1165 
1166  if ( eInd <= 3 ) {
1167  if ( elem2->giveNeighbor(1) == elem1->giveNumber() ) {
1168  ngb = eInd + 1;
1169  } else {
1170  ngb = 1;
1171  }
1172  } else {
1173  if ( elem2->giveNeighbor(eInd - 2) == elem1->giveNumber() ) {
1174  ngb = ( eInd > 4 ) ? eInd - 3 : 4;
1175  } else {
1176  ngb = eInd - 2;
1177  }
1178  }
1179 
1180  elem1 = elem2;
1181  }
1182 
1183  if ( elem2 != this ) {
1184 #ifdef DEBUG_CHECK
1185  #ifdef THREEPBT_3D
1186  if ( coords.at(1) > 0.000001 && coords.at(1) < 1999.99999 &&
1187  coords.at(2) > 0.000001 && coords.at(2) < 249.99999 &&
1188  coords.at(3) > 0.000001 && coords.at(3) < 499.99999 ) {
1189  if ( 987.5 - coords.at(1) > 0.000001 || coords.at(1) - 1012.5 > 0.000001 || 300.0 - coords.at(3) > 0.000001 ) {
1190  OOFEM_ERROR("Irregular %d [%d %d] not on boundary", iNum, iNode, jNode);
1191  }
1192  }
1193 
1194  #endif
1195 #endif
1196 #ifndef QUICK_HACK
1197  boundary = true;
1198 #endif
1199  // edge is on outer boundary
1200 
1201  // I do rely on the fact that if the list of connected elements is not availale
1202  // then the node is on top level
1203  iElems = mesh->giveNode(iNode)->giveConnectedElements();
1204  if ( !iElems->giveSize() ) {
1206  }
1207 
1208  jElems = mesh->giveNode(jNode)->giveConnectedElements();
1209  if ( !jElems->giveSize() ) {
1211  }
1212 
1213  IntArray common;
1214  if ( iElems->giveSize() <= jElems->giveSize() ) {
1215  common.preallocate( iElems->giveSize() );
1216  } else {
1217  common.preallocate( jElems->giveSize() );
1218  }
1219 
1220  // I do rely on the fact that the arrays are ordered !!!
1221  // I am using zero chunk because common is large enough
1222  elems = iElems->findCommonValuesSorted(* jElems, common, 0);
1223 #ifdef DEBUG_CHECK
1224  if ( !elems ) {
1225  OOFEM_ERROR("local outer edge %d %d is not shared by common elements",
1226  iNode, jNode);
1227  }
1228 
1229 #endif
1230  // after subdivision there will be twice as much of connected local elements
1231  irregular->preallocateConnectedElements(elems * 2);
1232  irregular->setNumber(-iNum); // mark local unshared irregular for connectivity setup
1233  // put the new node on appropriate edge of all local elements sharing both nodes
1234  // (if not yet done during neighbour traversal)
1235  for ( j = 1; j <= elems; j++ ) {
1236  elem = mesh->giveElement( common.at(j) );
1237  if ( elem == this ) {
1238  continue;
1239  }
1240 
1241  eInd = elem->giveEdgeIndex(iNode, jNode);
1242  if ( !elem->giveIrregular(eInd) ) {
1243  elem->setIrregular(eInd, iNum);
1244 
1245 #ifdef DEBUG_INFO
1246  #ifdef __PARALLEL_MODE
1247  // do not print global numbers of elements because they are not available (they are assigned at once after bisection);
1248  // do not print global numbers of irregulars as these may not be available yet
1249  OOFEM_LOG_INFO( "[%d] Irregular %d added on %d [%d] (edge %d, nodes %d %d [%d %d], nds %d %d %d %d [%d %d %d %d], ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
1250  mesh->giveSubdivision()->giveRank(), iNum,
1251  elem->giveNumber(), this->giveGlobalNumber(), eInd, iNode, jNode,
1253  elem->giveNode(1), elem->giveNode(2), elem->giveNode(3), elem->giveNode(4),
1254  mesh->giveNode( elem->giveNode(1) )->giveGlobalNumber(),
1255  mesh->giveNode( elem->giveNode(2) )->giveGlobalNumber(),
1256  mesh->giveNode( elem->giveNode(3) )->giveGlobalNumber(),
1257  mesh->giveNode( elem->giveNode(4) )->giveGlobalNumber(),
1258  elem->giveNeighbor(1), elem->giveNeighbor(2), elem->giveNeighbor(3), elem->giveNeighbor(4),
1259  elem->giveIrregular(1), elem->giveIrregular(2), elem->giveIrregular(3),
1260  elem->giveIrregular(4), elem->giveIrregular(5), elem->giveIrregular(6) );
1261  #else
1262  OOFEM_LOG_INFO( "Irregular %d added on %d (edge %d, nodes %d %d, nds %d %d %d %d, ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
1263  iNum, elem->giveNumber(), eInd, iNode, jNode,
1264  elem->giveNode(1), elem->giveNode(2), elem->giveNode(3), elem->giveNode(4),
1265  elem->giveNeighbor(1), elem->giveNeighbor(2), elem->giveNeighbor(3), elem->giveNeighbor(4),
1266  elem->giveIrregular(1), elem->giveIrregular(2), elem->giveIrregular(3),
1267  elem->giveIrregular(4), elem->giveIrregular(5), elem->giveIrregular(6) );
1268  #endif
1269 #endif
1270 
1271  if ( !elem->giveQueueFlag() ) {
1272  // add elem to list of elements for subdivision
1273  subdivqueue.push( elem->giveNumber() );
1274  elem->setQueueFlag(true);
1275  }
1276  }
1277  }
1278  }
1279  }
1280 
1281  if ( boundary ) {
1282  // OOFEM_LOG_INFO("Irregular %d set boundary\n", abs(irregular->giveNumber()));
1283  irregular->setBoundary(true);
1284  }
1285  }
1286  }
1287 
1288  this->setQueueFlag(false);
1289 }
1290 
1291 
1292 void
1294 {
1295 #ifdef DEBUG_CHECK
1296  if ( this->queue_flag ) {
1297  OOFEM_ERROR("unexpected queue flag of %d ", this->number);
1298  }
1299 
1300 #endif
1301 
1302  /* generates the children elements of already bisected element and adds them into mesh;
1303  * also children array of receiver is updated to contain children numbers */
1305  // no subdivision of receiver required
1306  children.clear();
1307  } else {
1308  int i, ind, nIrregulars = 0;
1309  int childNum, iedge, jedge, kedge, inode, jnode, knode;
1310  IntArray _nodes(3);
1313 #ifdef __PARALLEL_MODE
1314  bool ishared = false, jshared = false, kshared = false;
1315 #endif
1316 
1317  for ( i = 1; i <= irregular_nodes.giveSize(); i++ ) {
1318  if ( irregular_nodes.at(i) ) {
1319  nIrregulars++;
1320  }
1321  }
1322 
1323  this->children.resize(nIrregulars + 1);
1324 
1325  // leIndex determines primary edge to be subdivided
1326  /*
1327  * inode
1328  * o
1329  *
1330  *
1331  * kedge x x jedge
1332  *
1333  *
1334  * o x o
1335  * jnode iedge knode
1336  * =leIndex
1337  */
1338  iedge = leIndex;
1339  jedge = ( iedge < 3 ) ? iedge + 1 : 1;
1340  kedge = ( jedge < 3 ) ? jedge + 1 : 1;
1341 
1342  inode = ( iedge > 1 ) ? iedge - 1 : 3;
1343  jnode = ( inode < 3 ) ? inode + 1 : 1;
1344  knode = ( jnode < 3 ) ? jnode + 1 : 1;
1345 
1346 #ifdef __PARALLEL_MODE
1347  if ( shared_edges.giveSize() ) {
1348  if ( shared_edges.at(iedge) ) {
1349  ishared = true;
1350  }
1351 
1352  if ( shared_edges.at(jedge) ) {
1353  jshared = true;
1354  }
1355 
1356  if ( shared_edges.at(kedge) ) {
1357  kshared = true;
1358  }
1359  }
1360 
1361 #endif
1362 
1363  if ( irregular_nodes.at(iedge) && irregular_nodes.at(jedge) && irregular_nodes.at(kedge) ) {
1364  _nodes.at(1) = nodes.at(inode);
1365  _nodes.at(2) = irregular_nodes.at(kedge);
1366  _nodes.at(3) = irregular_nodes.at(iedge);
1367  childNum = mesh->giveNumberOfElements() + 1;
1368  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1369  mesh->addElement(child);
1370  children.at(1) = childNum;
1371  // set neigbour info
1372  child->setNeighbor( 1, -this->giveNeighbor(kedge) );
1373  child->setNeighbor(2, childNum + 2);
1374  child->setNeighbor(3, childNum + 1);
1375 #ifdef __PARALLEL_MODE
1376  // set shared info
1377  if ( kshared ) {
1378  child->makeSharedEdges();
1379  child->setSharedEdge( 1, shared_edges.at(kedge) );
1380  }
1381 
1382 #endif
1383 
1384  _nodes.at(1) = nodes.at(inode);
1385  _nodes.at(2) = irregular_nodes.at(iedge);
1386  _nodes.at(3) = irregular_nodes.at(jedge);
1387  childNum = mesh->giveNumberOfElements() + 1;
1388  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1389  mesh->addElement(child);
1390  children.at(2) = childNum;
1391  // set neigbour info
1392  child->setNeighbor(1, childNum - 1);
1393  child->setNeighbor(2, childNum + 2);
1394  child->setNeighbor( 3, -this->giveNeighbor(jedge) );
1395 #ifdef __PARALLEL_MODE
1396  // set shared info
1397  if ( jshared ) {
1398  child->makeSharedEdges();
1399  child->setSharedEdge( 3, shared_edges.at(jedge) );
1400  }
1401 
1402 #endif
1403 
1404  _nodes.at(1) = nodes.at(jnode);
1405  _nodes.at(2) = irregular_nodes.at(iedge);
1406  _nodes.at(3) = irregular_nodes.at(kedge);
1407  childNum = mesh->giveNumberOfElements() + 1;
1408  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1409  mesh->addElement(child);
1410  children.at(3) = childNum;
1411  // set neigbour info
1412  child->setNeighbor( 1, -this->giveNeighbor(iedge) );
1413  child->setNeighbor(2, childNum - 2);
1414  child->setNeighbor( 3, -this->giveNeighbor(kedge) );
1415 #ifdef __PARALLEL_MODE
1416  // set shared info
1417  if ( ishared || kshared ) {
1418  child->makeSharedEdges();
1419  if ( ishared ) {
1420  child->setSharedEdge( 1, shared_edges.at(iedge) );
1421  }
1422 
1423  if ( kshared ) {
1424  child->setSharedEdge( 3, shared_edges.at(kedge) );
1425  }
1426  }
1427 
1428 #endif
1429 
1430  _nodes.at(1) = nodes.at(knode);
1431  _nodes.at(2) = irregular_nodes.at(jedge);
1432  _nodes.at(3) = irregular_nodes.at(iedge);
1433  childNum = mesh->giveNumberOfElements() + 1;
1434  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1435  mesh->addElement(child);
1436  children.at(4) = childNum;
1437  // set neigbour info
1438  child->setNeighbor( 1, -this->giveNeighbor(jedge) );
1439  child->setNeighbor(2, childNum - 2);
1440  child->setNeighbor( 3, -this->giveNeighbor(iedge) );
1441 
1442 #ifdef __PARALLEL_MODE
1443  // set shared info
1444  if ( ishared || jshared ) {
1445  child->makeSharedEdges();
1446  if ( ishared ) {
1447  child->setSharedEdge( 3, shared_edges.at(iedge) );
1448  }
1449 
1450  if ( jshared ) {
1451  child->setSharedEdge( 1, shared_edges.at(jedge) );
1452  }
1453  }
1454 
1455 #endif
1456 
1457 #ifdef __PARALLEL_MODE
1458  // set connected elements
1459  if ( mesh->giveNode( nodes.at(inode) )->giveParallelMode() == DofManager_shared ) {
1460  mesh->giveNode( nodes.at(inode) )->eraseConnectedElement( this->giveNumber() );
1461  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 2);
1462  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 3);
1463  }
1464 
1465  if ( mesh->giveNode( nodes.at(jnode) )->giveParallelMode() == DofManager_shared ) {
1466  mesh->giveNode( nodes.at(jnode) )->eraseConnectedElement( this->giveNumber() );
1467  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum - 1);
1468  }
1469 
1470  if ( mesh->giveNode( nodes.at(knode) )->giveParallelMode() == DofManager_shared ) {
1471  mesh->giveNode( nodes.at(knode) )->eraseConnectedElement( this->giveNumber() );
1472  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum);
1473  }
1474 
1475  if ( ishared ) {
1476  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum);
1477  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 1);
1478  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 2);
1479  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 3);
1480  }
1481 
1482  if ( jshared ) {
1483  mesh->giveNode( irregular_nodes.at(jedge) )->insertConnectedElement(childNum);
1484  mesh->giveNode( irregular_nodes.at(jedge) )->insertConnectedElement(childNum - 2);
1485  }
1486 
1487  if ( kshared ) {
1488  mesh->giveNode( irregular_nodes.at(kedge) )->insertConnectedElement(childNum - 1);
1489  mesh->giveNode( irregular_nodes.at(kedge) )->insertConnectedElement(childNum - 3);
1490  }
1491 
1492 #endif
1493  } else if ( irregular_nodes.at(iedge) && irregular_nodes.at(jedge) ) {
1494  _nodes.at(1) = nodes.at(inode);
1495  _nodes.at(2) = nodes.at(jnode);
1496  _nodes.at(3) = irregular_nodes.at(iedge);
1497  childNum = mesh->giveNumberOfElements() + 1;
1498  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1499  mesh->addElement(child);
1500  children.at(1) = childNum;
1501  // set neigbour info
1502  child->setNeighbor( 1, -this->giveNeighbor(kedge) );
1503  child->setNeighbor( 2, -this->giveNeighbor(iedge) );
1504  child->setNeighbor(3, childNum + 1);
1505 #ifdef __PARALLEL_MODE
1506  // set shared info
1507  if ( ishared || kshared ) {
1508  child->makeSharedEdges();
1509  if ( ishared ) {
1510  child->setSharedEdge( 2, shared_edges.at(iedge) );
1511  }
1512 
1513  if ( kshared ) {
1514  child->setSharedEdge( 1, shared_edges.at(kedge) );
1515  }
1516  }
1517 
1518 #endif
1519 
1520  _nodes.at(1) = nodes.at(inode);
1521  _nodes.at(2) = irregular_nodes.at(iedge);
1522  _nodes.at(3) = irregular_nodes.at(jedge);
1523  childNum = mesh->giveNumberOfElements() + 1;
1524  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1525  mesh->addElement(child);
1526  children.at(2) = childNum;
1527  // set neigbour info
1528  child->setNeighbor(1, childNum - 1);
1529  child->setNeighbor(2, childNum + 1);
1530  child->setNeighbor( 3, -this->giveNeighbor(jedge) );
1531 #ifdef __PARALLEL_MODE
1532  // set shared info
1533  if ( jshared ) {
1534  child->makeSharedEdges();
1535  child->setSharedEdge( 3, shared_edges.at(jedge) );
1536  }
1537 
1538 #endif
1539 
1540  _nodes.at(1) = nodes.at(knode);
1541  _nodes.at(2) = irregular_nodes.at(jedge);
1542  _nodes.at(3) = irregular_nodes.at(iedge);
1543  childNum = mesh->giveNumberOfElements() + 1;
1544  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1545  mesh->addElement(child);
1546  children.at(3) = childNum;
1547  // set neigbour info
1548  child->setNeighbor( 1, -this->giveNeighbor(jedge) );
1549  child->setNeighbor(2, childNum - 1);
1550  child->setNeighbor( 3, -this->giveNeighbor(iedge) );
1551 #ifdef __PARALLEL_MODE
1552  // set shared info
1553  if ( ishared || jshared ) {
1554  child->makeSharedEdges();
1555  if ( ishared ) {
1556  child->setSharedEdge( 3, shared_edges.at(iedge) );
1557  }
1558 
1559  if ( jshared ) {
1560  child->setSharedEdge( 1, shared_edges.at(jedge) );
1561  }
1562  }
1563 
1564 #endif
1565 
1566 #ifdef __PARALLEL_MODE
1567  // set connected elements
1568  if ( mesh->giveNode( nodes.at(inode) )->giveParallelMode() == DofManager_shared ) {
1569  mesh->giveNode( nodes.at(inode) )->eraseConnectedElement( this->giveNumber() );
1570  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 1);
1571  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 2);
1572  }
1573 
1574  if ( mesh->giveNode( nodes.at(jnode) )->giveParallelMode() == DofManager_shared ) {
1575  mesh->giveNode( nodes.at(jnode) )->eraseConnectedElement( this->giveNumber() );
1576  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum - 2);
1577  }
1578 
1579  if ( mesh->giveNode( nodes.at(knode) )->giveParallelMode() == DofManager_shared ) {
1580  mesh->giveNode( nodes.at(knode) )->eraseConnectedElement( this->giveNumber() );
1581  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum);
1582  }
1583 
1584  if ( ishared ) {
1585  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum);
1586  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 1);
1587  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 2);
1588  }
1589 
1590  if ( jshared ) {
1591  mesh->giveNode( irregular_nodes.at(jedge) )->insertConnectedElement(childNum);
1592  mesh->giveNode( irregular_nodes.at(jedge) )->insertConnectedElement(childNum - 1);
1593  }
1594 
1595 #endif
1596  } else if ( irregular_nodes.at(iedge) && irregular_nodes.at(kedge) ) {
1597  _nodes.at(1) = nodes.at(inode);
1598  _nodes.at(2) = irregular_nodes.at(kedge);
1599  _nodes.at(3) = irregular_nodes.at(iedge);
1600  childNum = mesh->giveNumberOfElements() + 1;
1601  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1602  mesh->addElement(child);
1603  children.at(1) = childNum;
1604  // set neigbour info
1605  child->setNeighbor( 1, -this->giveNeighbor(kedge) );
1606  child->setNeighbor(2, childNum + 1);
1607  child->setNeighbor(3, childNum + 2);
1608 #ifdef __PARALLEL_MODE
1609  // set shared info
1610  if ( kshared ) {
1611  child->makeSharedEdges();
1612  child->setSharedEdge( 1, shared_edges.at(kedge) );
1613  }
1614 
1615 #endif
1616 
1617  _nodes.at(1) = nodes.at(jnode);
1618  _nodes.at(2) = irregular_nodes.at(iedge);
1619  _nodes.at(3) = irregular_nodes.at(kedge);
1620  childNum = mesh->giveNumberOfElements() + 1;
1621  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1622  mesh->addElement(child);
1623  children.at(2) = childNum;
1624  // set neigbour info
1625  child->setNeighbor( 1, -this->giveNeighbor(iedge) );
1626  child->setNeighbor(2, childNum - 1);
1627  child->setNeighbor( 3, -this->giveNeighbor(kedge) );
1628 #ifdef __PARALLEL_MODE
1629  // set shared info
1630  if ( ishared || kshared ) {
1631  child->makeSharedEdges();
1632  if ( ishared ) {
1633  child->setSharedEdge( 1, shared_edges.at(iedge) );
1634  }
1635 
1636  if ( kshared ) {
1637  child->setSharedEdge( 3, shared_edges.at(kedge) );
1638  }
1639  }
1640 
1641 #endif
1642 
1643  _nodes.at(1) = nodes.at(inode);
1644  _nodes.at(2) = irregular_nodes.at(iedge);
1645  _nodes.at(3) = nodes.at(knode);
1646  childNum = mesh->giveNumberOfElements() + 1;
1647  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1648  mesh->addElement(child);
1649  children.at(3) = childNum;
1650  // set neigbour info
1651  child->setNeighbor(1, childNum - 2);
1652  child->setNeighbor( 2, -this->giveNeighbor(iedge) );
1653  child->setNeighbor( 3, -this->giveNeighbor(jedge) );
1654 #ifdef __PARALLEL_MODE
1655  // set shared info
1656  if ( ishared || jshared ) {
1657  child->makeSharedEdges();
1658  if ( ishared ) {
1659  child->setSharedEdge( 2, shared_edges.at(iedge) );
1660  }
1661 
1662  if ( jshared ) {
1663  child->setSharedEdge( 3, shared_edges.at(jedge) );
1664  }
1665  }
1666 
1667 #endif
1668 
1669 #ifdef __PARALLEL_MODE
1670  // set connected elements
1671  if ( mesh->giveNode( nodes.at(inode) )->giveParallelMode() == DofManager_shared ) {
1672  mesh->giveNode( nodes.at(inode) )->eraseConnectedElement( this->giveNumber() );
1673  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum);
1674  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 2);
1675  }
1676 
1677  if ( mesh->giveNode( nodes.at(jnode) )->giveParallelMode() == DofManager_shared ) {
1678  mesh->giveNode( nodes.at(jnode) )->eraseConnectedElement( this->giveNumber() );
1679  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum - 1);
1680  }
1681 
1682  if ( mesh->giveNode( nodes.at(knode) )->giveParallelMode() == DofManager_shared ) {
1683  mesh->giveNode( nodes.at(knode) )->eraseConnectedElement( this->giveNumber() );
1684  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum);
1685  }
1686 
1687  if ( ishared ) {
1688  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum);
1689  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 1);
1690  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 2);
1691  }
1692 
1693  if ( kshared ) {
1694  mesh->giveNode( irregular_nodes.at(kedge) )->insertConnectedElement(childNum - 1);
1695  mesh->giveNode( irregular_nodes.at(kedge) )->insertConnectedElement(childNum - 2);
1696  }
1697 
1698 #endif
1699  } else if ( irregular_nodes.at(iedge) ) {
1700  _nodes.at(1) = nodes.at(inode);
1701  _nodes.at(2) = nodes.at(jnode);
1702  _nodes.at(3) = irregular_nodes.at(iedge);
1703  childNum = mesh->giveNumberOfElements() + 1;
1704  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1705  mesh->addElement(child);
1706  children.at(1) = childNum;
1707  // set neigbour info
1708  child->setNeighbor( 1, -this->giveNeighbor(kedge) );
1709  child->setNeighbor( 2, -this->giveNeighbor(iedge) );
1710  child->setNeighbor(3, childNum + 1);
1711 #ifdef __PARALLEL_MODE
1712  // set shared info
1713  if ( ishared || kshared ) {
1714  child->makeSharedEdges();
1715  if ( ishared ) {
1716  child->setSharedEdge( 2, shared_edges.at(iedge) );
1717  }
1718 
1719  if ( kshared ) {
1720  child->setSharedEdge( 1, shared_edges.at(kedge) );
1721  }
1722  }
1723 
1724 #endif
1725 
1726  _nodes.at(1) = nodes.at(inode);
1727  _nodes.at(2) = irregular_nodes.at(iedge);
1728  _nodes.at(3) = nodes.at(knode);
1729  childNum = mesh->giveNumberOfElements() + 1;
1730  child = new Subdivision :: RS_Triangle(childNum, mesh, this->number, _nodes); // number, parent, coords
1731  mesh->addElement(child);
1732  children.at(2) = childNum;
1733  // set neigbour info
1734  child->setNeighbor(1, childNum - 1);
1735  child->setNeighbor( 2, -this->giveNeighbor(iedge) );
1736  child->setNeighbor( 3, -this->giveNeighbor(jedge) );
1737 #ifdef __PARALLEL_MODE
1738  // set shared info
1739  if ( ishared || jshared ) {
1740  child->makeSharedEdges();
1741  if ( ishared ) {
1742  child->setSharedEdge( 2, shared_edges.at(iedge) );
1743  }
1744 
1745  if ( jshared ) {
1746  child->setSharedEdge( 3, shared_edges.at(jedge) );
1747  }
1748  }
1749 
1750 #endif
1751 
1752 #ifdef __PARALLEL_MODE
1753  // set connected elements
1754  if ( mesh->giveNode( nodes.at(inode) )->giveParallelMode() == DofManager_shared ) {
1755  mesh->giveNode( nodes.at(inode) )->eraseConnectedElement( this->giveNumber() );
1756  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum);
1757  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 1);
1758  }
1759 
1760  if ( mesh->giveNode( nodes.at(jnode) )->giveParallelMode() == DofManager_shared ) {
1761  mesh->giveNode( nodes.at(jnode) )->eraseConnectedElement( this->giveNumber() );
1762  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum - 1);
1763  }
1764 
1765  if ( mesh->giveNode( nodes.at(knode) )->giveParallelMode() == DofManager_shared ) {
1766  mesh->giveNode( nodes.at(knode) )->eraseConnectedElement( this->giveNumber() );
1767  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum);
1768  }
1769 
1770  if ( ishared ) {
1771  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum);
1772  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 1);
1773  }
1774 
1775 #endif
1776  } else {
1777  OOFEM_ERROR("element %d internal data inconsistency", this->number);
1778  }
1779 
1780  // if there is neighbor of "this" not designated for bisection
1781  // its neihgbor corresponding to "this" must be made negative to enforce update_neighbors
1782  for ( i = 1; i <= 3; i++ ) {
1783  if ( this->neghbours_base_elements.at(i) ) {
1784 #ifdef DEBUG_CHECK
1785  if ( this->neghbours_base_elements.at(i) < 0 ) {
1786  OOFEM_ERROR("negative neighbor of %d not expected", this->number);
1787  }
1788 
1789 #endif
1790 
1791  ngb = mesh->giveElement( this->neghbours_base_elements.at(i) );
1792  if ( !ngb->hasIrregulars() ) {
1793  ind = ngb->giveNeighbors()->findFirstIndexOf(this->number);
1794  if ( ind ) {
1795  ngb->setNeighbor(ind, -this->number);
1796  }
1797  }
1798  }
1799  }
1800  }
1801 }
1802 
1803 
1804 void
1806 {
1807 #ifdef DEBUG_CHECK
1808  if ( this->queue_flag ) {
1809  OOFEM_ERROR("unexpected queue flag of %d ", this->number);
1810  }
1811 
1812 #endif
1813 
1814  /* generates the children elements of already bisected element and adds them into mesh;
1815  * this is done recursively;
1816  * also children array of receiver is updated to contain children numbers */
1817 
1819  // no subdivision of receiver required
1820  children.clear();
1821  } else {
1822  int irregulars1 = 0, irregulars2 = 0;
1823  int childNum, iedge, jedge, kedge, iiedge, jjedge, kkedge, inode, jnode, knode, nnode, iside, jside, kside, nside;
1824  int i, ind, leIndex1 = 0, leIndex2 = 0;
1825  IntArray _nodes(4);
1826  Subdivision :: RS_Tetra *child;
1828 #ifdef __PARALLEL_MODE
1829  int eNum = this->mesh->giveNumberOfEdges();
1831  bool ishared = false, jshared = false, kshared = false;
1832  bool iishared = false, jjshared = false, kkshared = false;
1833 #endif
1834 
1835  // leIndex determines primary edge to be subdivided
1836  // it is identified either with iedge or iiedge
1837  /*
1838  *
1839  * inode
1840  * o
1841  *
1842  * iiedge=leIndex(if>3)
1843  * x
1844  * inode, jnode, knode - base nodes
1845  * kedge x nnode x jedge nnode - top node
1846  * o iiedge, jjedge, kkedge - edges connecting base nodes with top node
1847  *
1848  * jjedge x x kkedge
1849  *
1850  * o x o
1851  * jnode iedge knode
1852  * =leIndex(if<=3)
1853  *
1854  * children keep all nodes in the same order except for that one which is replaced by irregular on leIndex edge !!!
1855  */
1856  if ( leIndex <= 3 ) {
1857  iedge = leIndex;
1858  } else {
1859  iedge = ( leIndex == 6 ) ? 1 : leIndex - 2;
1860  }
1861 
1862  jedge = ( iedge < 3 ) ? iedge + 1 : 1;
1863  kedge = ( jedge < 3 ) ? jedge + 1 : 1;
1864 
1865  inode = ( iedge > 1 ) ? iedge - 1 : 3;
1866  jnode = ( inode < 3 ) ? inode + 1 : 1;
1867  knode = ( jnode < 3 ) ? jnode + 1 : 1;
1868  nnode = 4;
1869 
1870  iiedge = inode + 3;
1871  jjedge = jnode + 3;
1872  kkedge = knode + 3;
1873 
1874  iside = iedge + 1;
1875  jside = jedge + 1;
1876  kside = kedge + 1;
1877  nside = 1;
1878 
1879  this->children.resize(2);
1880 
1881 #ifdef __PARALLEL_MODE
1882  if ( shared_edges.giveSize() ) {
1883  if ( shared_edges.at(iedge) ) {
1884  ishared = true;
1885  }
1886 
1887  if ( shared_edges.at(jedge) ) {
1888  jshared = true;
1889  }
1890 
1891  if ( shared_edges.at(kedge) ) {
1892  kshared = true;
1893  }
1894 
1895  if ( shared_edges.at(iiedge) ) {
1896  iishared = true;
1897  }
1898 
1899  if ( shared_edges.at(jjedge) ) {
1900  jjshared = true;
1901  }
1902 
1903  if ( shared_edges.at(kkedge) ) {
1904  kkshared = true;
1905  }
1906  }
1907 
1908 #endif
1909 
1910  if ( leIndex <= 3 ) {
1911  // count number of irregulars on each child;
1912  // if there is one irregular only, the set leIndex is valid
1913  // otherwise it is corrected later
1914  if ( irregular_nodes.at(iiedge) ) {
1915  irregulars1++;
1916  irregulars2++;
1917  leIndex1 = leIndex2 = 4;
1918  }
1919 
1920  if ( irregular_nodes.at(kedge) ) {
1921  irregulars1++;
1922  leIndex1 = 1;
1923  }
1924 
1925  if ( irregular_nodes.at(jjedge) ) {
1926  irregulars1++;
1927  leIndex1 = 5;
1928  }
1929 
1930  if ( irregular_nodes.at(jedge) ) {
1931  irregulars2++;
1932  leIndex2 = 3;
1933  }
1934 
1935  if ( irregular_nodes.at(kkedge) ) {
1936  irregulars2++;
1937  leIndex2 = 6;
1938  }
1939 
1940  if ( irregulars1 > 1 ) {
1941  // use parent information about the longest edge on kside
1942  leIndex1 = this->side_leIndex.at(kside);
1943 
1944  if ( leIndex1 == jjedge ) {
1945  leIndex1 = 5;
1946  } else if ( leIndex1 == kedge ) {
1947  leIndex1 = 1;
1948  } else {
1949 #ifdef DEBUG_CHECK
1950  if ( leIndex1 != iiedge ) {
1951  OOFEM_ERROR("side longest edge inconsistency on %d", this->number);
1952  }
1953 
1954 #endif
1955  leIndex1 = 4;
1956  }
1957  }
1958 
1959  if ( irregulars2 > 1 ) {
1960  // use parent information about the longest edge on jside
1961  leIndex2 = this->side_leIndex.at(jside);
1962 
1963  if ( leIndex2 == kkedge ) {
1964  leIndex2 = 6;
1965  } else if ( leIndex2 == jedge ) {
1966  leIndex2 = 3;
1967  } else {
1968 #ifdef DEBUG_CHECK
1969  if ( leIndex2 != iiedge ) {
1970  OOFEM_ERROR("side longest edge inconsistency on %d", this->number);
1971  }
1972 
1973 #endif
1974  leIndex2 = 4;
1975  }
1976  }
1977 
1978 #ifdef __PARALLEL_MODE
1979  int i_shared_id = 0, n_shared_id = 0;
1980 
1981  // check whether new edges are potentially shared
1982  if ( ishared ) {
1983  if ( mesh->giveNode( nodes.at(inode) )->giveParallelMode() == DofManager_shared ) {
1984  if ( !this->giveNeighbor(nside) ) {
1985  _edge = new Subdivision :: RS_SharedEdge(this->mesh);
1986  _edge->setEdgeNodes( irregular_nodes.at(iedge), nodes.at(inode) );
1987 
1988  this->mesh->addEdge(_edge);
1989  eNum++;
1990 
1991  // get the tentative partitions
1992  // they must be confirmed via mutual communication
1993  IntArray partitions;
1994  if ( _edge->giveSharedPartitions(partitions) ) {
1995  i_shared_id = eNum;
1996  _edge->setPartitions(partitions);
1997  // put edge number into queue of shared edges to resolve remote partitions
1998  sharedEdgesQueue.push_back(eNum);
1999  }
2000  }
2001  }
2002 
2003  if ( mesh->giveNode( nodes.at(nnode) )->giveParallelMode() == DofManager_shared ) {
2004  if ( !this->giveNeighbor(iside) ) {
2005  _edge = new Subdivision :: RS_SharedEdge(this->mesh);
2006  _edge->setEdgeNodes( irregular_nodes.at(iedge), nodes.at(nnode) );
2007 
2008  this->mesh->addEdge(_edge);
2009  eNum++;
2010 
2011  // get the tentative partitions
2012  // they must be confirmed via mutual communication
2013  IntArray partitions;
2014  if ( _edge->giveSharedPartitions(partitions) ) {
2015  n_shared_id = eNum;
2016  _edge->setPartitions(partitions);
2017  // put edge number into queue of shared edges to resolve remote partitions
2018  sharedEdgesQueue.push_back(eNum);
2019  }
2020  }
2021  }
2022  }
2023 
2024 #endif
2025 
2026  _nodes.at(1) = nodes.at(inode);
2027  _nodes.at(2) = nodes.at(jnode);
2028  _nodes.at(3) = irregular_nodes.at(iedge);
2029  _nodes.at(4) = nodes.at(nnode);
2030  childNum = mesh->giveNumberOfElements() + 1;
2031  child = new Subdivision :: RS_Tetra(childNum, mesh, this->number, _nodes); // number, parent, coords
2032  mesh->addElement(child);
2033  children.at(1) = childNum;
2034 
2035  // set neigbour info
2036  if ( irregulars1 ) {
2037  child->setNeighbor( 1, this->giveNeighbor(nside) );
2038  child->setNeighbor( 2, this->giveNeighbor(kside) );
2039  child->setNeighbor( 3, this->giveNeighbor(iside) );
2040 
2041  child->setIrregular( 1, irregular_nodes.at(kedge) );
2042  child->setIrregular( 5, irregular_nodes.at(jjedge) );
2043  child->setIrregular( 4, irregular_nodes.at(iiedge) );
2044  } else {
2045  child->setNeighbor( 1, -this->giveNeighbor(nside) );
2046  child->setNeighbor( 2, -this->giveNeighbor(kside) );
2047  child->setNeighbor( 3, -this->giveNeighbor(iside) );
2048  }
2049 
2050  // neihgbor4 of child1 is changed to negative during subdivision (if any) of child2
2051  child->setNeighbor(4, childNum + 1);
2052 
2053 #ifdef __PARALLEL_MODE
2054  // set shared info
2055  if ( ishared || kshared || iishared || jjshared ) {
2056  child->makeSharedEdges();
2057  if ( ishared ) {
2058  child->setSharedEdge( 2, shared_edges.at(iedge) );
2059  }
2060 
2061  if ( kshared ) {
2062  child->setSharedEdge( 1, shared_edges.at(kedge) );
2063  }
2064 
2065  if ( iishared ) {
2066  child->setSharedEdge( 4, shared_edges.at(iiedge) );
2067  }
2068 
2069  if ( jjshared ) {
2070  child->setSharedEdge( 5, shared_edges.at(jjedge) );
2071  }
2072 
2073  child->setSharedEdge(3, i_shared_id);
2074  child->setSharedEdge(6, n_shared_id);
2075  }
2076 
2077 #endif
2078 
2079 #ifdef DEBUG_INFO
2080  #ifdef __PARALLEL_MODE
2081  OOFEM_LOG_INFO("[%d] Child %d generated on parent %d [%d] (leIndex %d, nds %d %d %d %d [%d %d %d %d], ngbs %d %d %d %d, irr %d %d %d %d %d %d [%d %d %d %d %d %d])\n",
2082  mesh->giveSubdivision()->giveRank(), childNum,
2083  this->number, this->giveGlobalNumber(), this->leIndex,
2084  _nodes.at(1), _nodes.at(2), _nodes.at(3), _nodes.at(4),
2085  mesh->giveNode( _nodes.at(1) )->giveGlobalNumber(),
2086  mesh->giveNode( _nodes.at(2) )->giveGlobalNumber(),
2087  mesh->giveNode( _nodes.at(3) )->giveGlobalNumber(),
2088  mesh->giveNode( _nodes.at(4) )->giveGlobalNumber(),
2089  child->giveNeighbor(1), child->giveNeighbor(2), child->giveNeighbor(3), child->giveNeighbor(4),
2090  child->giveIrregular(1), child->giveIrregular(2), child->giveIrregular(3),
2091  child->giveIrregular(4), child->giveIrregular(5), child->giveIrregular(6),
2092  ( child->giveIrregular(1) ) ? mesh->giveNode( child->giveIrregular(1) )->giveGlobalNumber() : 0,
2093  ( child->giveIrregular(2) ) ? mesh->giveNode( child->giveIrregular(2) )->giveGlobalNumber() : 0,
2094  ( child->giveIrregular(3) ) ? mesh->giveNode( child->giveIrregular(3) )->giveGlobalNumber() : 0,
2095  ( child->giveIrregular(4) ) ? mesh->giveNode( child->giveIrregular(4) )->giveGlobalNumber() : 0,
2096  ( child->giveIrregular(5) ) ? mesh->giveNode( child->giveIrregular(5) )->giveGlobalNumber() : 0,
2097  ( child->giveIrregular(6) ) ? mesh->giveNode( child->giveIrregular(6) )->giveGlobalNumber() : 0);
2098  #else
2099  OOFEM_LOG_INFO( "Child %d generated on parent %d (leIndex %d, nds %d %d %d %d, ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
2100  childNum, this->number, this->leIndex, _nodes.at(1), _nodes.at(2), _nodes.at(3), _nodes.at(4),
2101  child->giveNeighbor(1), child->giveNeighbor(2), child->giveNeighbor(3), child->giveNeighbor(4),
2102  child->giveIrregular(1), child->giveIrregular(2), child->giveIrregular(3),
2103  child->giveIrregular(4), child->giveIrregular(5), child->giveIrregular(6) );
2104  #endif
2105 #endif
2106 
2107  _nodes.at(1) = nodes.at(inode);
2108  _nodes.at(2) = irregular_nodes.at(iedge);
2109  _nodes.at(3) = nodes.at(knode);
2110  _nodes.at(4) = nodes.at(nnode);
2111  childNum = mesh->giveNumberOfElements() + 1;
2112  child = new Subdivision :: RS_Tetra(childNum, mesh, this->number, _nodes); // number, parent, coords
2113  mesh->addElement(child);
2114  children.at(2) = childNum;
2115 
2116  // set neigbour info
2117  if ( irregulars2 ) {
2118  child->setNeighbor( 1, this->giveNeighbor(nside) );
2119  child->setNeighbor( 3, this->giveNeighbor(iside) );
2120  child->setNeighbor( 4, this->giveNeighbor(jside) );
2121 
2122  child->setIrregular( 3, irregular_nodes.at(jedge) );
2123  child->setIrregular( 6, irregular_nodes.at(kkedge) );
2124  child->setIrregular( 4, irregular_nodes.at(iiedge) );
2125  } else {
2126  child->setNeighbor( 1, -this->giveNeighbor(nside) );
2127  child->setNeighbor( 3, -this->giveNeighbor(iside) );
2128  child->setNeighbor( 4, -this->giveNeighbor(jside) );
2129  }
2130 
2131  // neihgbor2 of child2 is changed to negative during subdivision (if any) of child1
2132  child->setNeighbor(2, childNum - 1);
2133 #ifdef __PARALLEL_MODE
2134  // set shared info
2135  if ( ishared || jshared || iishared || kkshared ) {
2136  child->makeSharedEdges();
2137  if ( ishared ) {
2138  child->setSharedEdge( 2, shared_edges.at(iedge) );
2139  }
2140 
2141  if ( jshared ) {
2142  child->setSharedEdge( 3, shared_edges.at(jedge) );
2143  }
2144 
2145  if ( iishared ) {
2146  child->setSharedEdge( 4, shared_edges.at(iiedge) );
2147  }
2148 
2149  if ( kkshared ) {
2150  child->setSharedEdge( 6, shared_edges.at(kkedge) );
2151  }
2152 
2153  child->setSharedEdge(1, i_shared_id);
2154  child->setSharedEdge(5, n_shared_id);
2155  }
2156 
2157 #endif
2158 
2159 #ifdef __PARALLEL_MODE
2160  // set connected elements
2161  if ( mesh->giveNode( nodes.at(inode) )->giveParallelMode() == DofManager_shared ) {
2162  mesh->giveNode( nodes.at(inode) )->eraseConnectedElement( this->giveNumber() );
2163  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum);
2164  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 1);
2165  }
2166 
2167  if ( mesh->giveNode( nodes.at(nnode) )->giveParallelMode() == DofManager_shared ) {
2168  mesh->giveNode( nodes.at(nnode) )->eraseConnectedElement( this->giveNumber() );
2169  mesh->giveNode( nodes.at(nnode) )->insertConnectedElement(childNum);
2170  mesh->giveNode( nodes.at(nnode) )->insertConnectedElement(childNum - 1);
2171  }
2172 
2173  if ( mesh->giveNode( nodes.at(jnode) )->giveParallelMode() == DofManager_shared ) {
2174  mesh->giveNode( nodes.at(jnode) )->eraseConnectedElement( this->giveNumber() );
2175  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum - 1);
2176  }
2177 
2178  if ( mesh->giveNode( nodes.at(knode) )->giveParallelMode() == DofManager_shared ) {
2179  mesh->giveNode( nodes.at(knode) )->eraseConnectedElement( this->giveNumber() );
2180  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum);
2181  }
2182 
2183  if ( ishared ) {
2184  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum);
2185  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 1);
2186  }
2187 
2188 #endif
2189 
2190 #ifdef DEBUG_INFO
2191  #ifdef __PARALLEL_MODE
2192  OOFEM_LOG_INFO("[%d] Child %d generated on parent %d [%d] (leIndex %d, nds %d %d %d %d [%d %d %d %d], ngbs %d %d %d %d, irr %d %d %d %d %d %d [%d %d %d %d %d %d])\n",
2193  mesh->giveSubdivision()->giveRank(), childNum,
2194  this->number, this->giveGlobalNumber(), this->leIndex,
2195  _nodes.at(1), _nodes.at(2), _nodes.at(3), _nodes.at(4),
2196  mesh->giveNode( _nodes.at(1) )->giveGlobalNumber(),
2197  mesh->giveNode( _nodes.at(2) )->giveGlobalNumber(),
2198  mesh->giveNode( _nodes.at(3) )->giveGlobalNumber(),
2199  mesh->giveNode( _nodes.at(4) )->giveGlobalNumber(),
2200  child->giveNeighbor(1), child->giveNeighbor(2), child->giveNeighbor(3), child->giveNeighbor(4),
2201  child->giveIrregular(1), child->giveIrregular(2), child->giveIrregular(3),
2202  child->giveIrregular(4), child->giveIrregular(5), child->giveIrregular(6),
2203  ( child->giveIrregular(1) ) ? mesh->giveNode( child->giveIrregular(1) )->giveGlobalNumber() : 0,
2204  ( child->giveIrregular(2) ) ? mesh->giveNode( child->giveIrregular(2) )->giveGlobalNumber() : 0,
2205  ( child->giveIrregular(3) ) ? mesh->giveNode( child->giveIrregular(3) )->giveGlobalNumber() : 0,
2206  ( child->giveIrregular(4) ) ? mesh->giveNode( child->giveIrregular(4) )->giveGlobalNumber() : 0,
2207  ( child->giveIrregular(5) ) ? mesh->giveNode( child->giveIrregular(5) )->giveGlobalNumber() : 0,
2208  ( child->giveIrregular(6) ) ? mesh->giveNode( child->giveIrregular(6) )->giveGlobalNumber() : 0);
2209  #else
2210  OOFEM_LOG_INFO( "Child %d generated on parent %d (leIndex %d, nds %d %d %d %d, ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
2211  childNum, this->number, this->leIndex, _nodes.at(1), _nodes.at(2), _nodes.at(3), _nodes.at(4),
2212  child->giveNeighbor(1), child->giveNeighbor(2), child->giveNeighbor(3), child->giveNeighbor(4),
2213  child->giveIrregular(1), child->giveIrregular(2), child->giveIrregular(3),
2214  child->giveIrregular(4), child->giveIrregular(5), child->giveIrregular(6) );
2215  #endif
2216 #endif
2217 
2218  // set connected elements to nonshared nodes on outer boundary
2219 #ifdef __PARALLEL_MODE
2220  if ( !ishared ) {
2221 #endif
2222  if ( mesh->giveNode( irregular_nodes.at(iedge) )->giveNumber() < 0 ) { // check for marked local irregular
2223  // irregular node on outer boundary
2224  // insert connected elements
2225  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum);
2226  mesh->giveNode( irregular_nodes.at(iedge) )->insertConnectedElement(childNum - 1);
2227  }
2228 
2229  // update connectivity of both end nodes of iedge if not shared
2230  // and ONLY if there is already list of connected elements
2231 #ifdef __PARALLEL_MODE
2232  if ( mesh->giveNode( nodes.at(jnode) )->giveParallelMode() != DofManager_shared ) { // update already done
2233 #endif
2234  if ( mesh->giveNode( nodes.at(jnode) )->giveConnectedElements()->giveSize() ) {
2235  mesh->giveNode( nodes.at(jnode) )->eraseConnectedElement( this->giveNumber() );
2236  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum - 1);
2237  }
2238 
2239 #ifdef __PARALLEL_MODE
2240  }
2241 
2242 #endif
2243 #ifdef __PARALLEL_MODE
2244  if ( mesh->giveNode( nodes.at(knode) )->giveParallelMode() != DofManager_shared ) { // update already done
2245 #endif
2246  if ( mesh->giveNode( nodes.at(knode) )->giveConnectedElements()->giveSize() ) {
2247  mesh->giveNode( nodes.at(knode) )->eraseConnectedElement( this->giveNumber() );
2248  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum);
2249  }
2250 
2251 #ifdef __PARALLEL_MODE
2252  }
2253 #endif
2254 #ifdef __PARALLEL_MODE
2255  }
2256 #endif
2257 #ifdef __PARALLEL_MODE
2258  if ( !iishared ) {
2259 #endif
2260  // update connectivity of both end nodes of iiedge if not shared
2261  // and ONLY if there is already list of connected elements
2262 #ifdef __PARALLEL_MODE
2263  if ( mesh->giveNode( nodes.at(inode) )->giveParallelMode() != DofManager_shared ) { // update already done
2264 #endif
2265  if ( mesh->giveNode( nodes.at(inode) )->giveConnectedElements()->giveSize() ) {
2266  mesh->giveNode( nodes.at(inode) )->eraseConnectedElement( this->giveNumber() );
2267  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum);
2268  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 1);
2269  }
2270 
2271 #ifdef __PARALLEL_MODE
2272  }
2273 #endif
2274 #ifdef __PARALLEL_MODE
2275  if ( mesh->giveNode( nodes.at(nnode) )->giveParallelMode() != DofManager_shared ) { // update already done
2276 #endif
2277  if ( mesh->giveNode( nodes.at(nnode) )->giveConnectedElements()->giveSize() ) {
2278  mesh->giveNode( nodes.at(nnode) )->eraseConnectedElement( this->giveNumber() );
2279  mesh->giveNode( nodes.at(nnode) )->insertConnectedElement(childNum);
2280  mesh->giveNode( nodes.at(nnode) )->insertConnectedElement(childNum - 1);
2281  }
2282 
2283 #ifdef __PARALLEL_MODE
2284  }
2285 #endif
2286 #ifdef __PARALLEL_MODE
2287  }
2288 #endif
2289  } else {
2290  // count number of irregulars on each child;
2291  // if there is one irregular only, the set leIndex is valid
2292  // otherwise it is corrected later
2293  if ( irregular_nodes.at(iedge) ) {
2294  irregulars1++;
2295  irregulars2++;
2296  leIndex1 = leIndex2 = 2;
2297  }
2298 
2299  if ( irregular_nodes.at(jedge) ) {
2300  irregulars1++;
2301  leIndex1 = 3;
2302  }
2303 
2304  if ( irregular_nodes.at(kedge) ) {
2305  irregulars1++;
2306  leIndex1 = 1;
2307  }
2308 
2309  if ( irregular_nodes.at(jjedge) ) {
2310  irregulars2++;
2311  leIndex2 = 5;
2312  }
2313 
2314  if ( irregular_nodes.at(kkedge) ) {
2315  irregulars2++;
2316  leIndex2 = 6;
2317  }
2318 
2319  if ( irregulars1 > 1 ) {
2320  // use parent information about the longest edge on nside
2321  leIndex1 = this->side_leIndex.at(nside);
2322 
2323  if ( leIndex1 == jedge ) {
2324  leIndex1 = 3;
2325  } else if ( leIndex1 == kedge ) {
2326  leIndex1 = 1;
2327  } else {
2328 #ifdef DEBUG_CHECK
2329  if ( leIndex1 != iedge ) {
2330  OOFEM_ERROR("side longest edge inconsistency on %d", this->number);
2331  }
2332 
2333 #endif
2334  leIndex1 = 2;
2335  }
2336  }
2337 
2338  if ( irregulars2 > 1 ) {
2339  // use parent information about the longest edge on iside
2340  leIndex2 = this->side_leIndex.at(iside);
2341 
2342  if ( leIndex2 == kkedge ) {
2343  leIndex2 = 6;
2344  } else if ( leIndex2 == jjedge ) {
2345  leIndex2 = 5;
2346  } else {
2347 #ifdef DEBUG_CHECK
2348  if ( leIndex2 != iedge ) {
2349  OOFEM_ERROR("side longest edge inconsistency on %d", this->number);
2350  }
2351 
2352 #endif
2353  leIndex2 = 2;
2354  }
2355  }
2356 
2357 #ifdef __PARALLEL_MODE
2358  int j_shared_id = 0, k_shared_id = 0;
2359 
2360  // check whether new edges are potentially shared
2361  if ( iishared ) {
2362  if ( mesh->giveNode( nodes.at(jnode) )->giveParallelMode() == DofManager_shared ) {
2363  if ( !this->giveNeighbor(kside) ) {
2364  _edge = new Subdivision :: RS_SharedEdge(this->mesh);
2365  _edge->setEdgeNodes( irregular_nodes.at(iiedge), nodes.at(jnode) );
2366 
2367  this->mesh->addEdge(_edge);
2368  eNum++;
2369 
2370  // get the tentative partitions
2371  // they must be confirmed via mutual communication
2372  IntArray partitions;
2373  if ( _edge->giveSharedPartitions(partitions) ) {
2374  j_shared_id = eNum;
2375  _edge->setPartitions(partitions);
2376  // put edge number into queue of shared edges to resolve remote partitions
2377  sharedEdgesQueue.push_back(eNum);
2378  }
2379  }
2380  }
2381 
2382  if ( mesh->giveNode( nodes.at(knode) )->giveParallelMode() == DofManager_shared ) {
2383  if ( !this->giveNeighbor(jside) ) {
2384  _edge = new Subdivision :: RS_SharedEdge(this->mesh);
2385  _edge->setEdgeNodes( irregular_nodes.at(iiedge), nodes.at(knode) );
2386 
2387  this->mesh->addEdge(_edge);
2388  eNum++;
2389 
2390  // get the tentative partitions
2391  // they must be confirmed via mutual communication
2392  IntArray partitions;
2393  if ( _edge->giveSharedPartitions(partitions) ) {
2394  k_shared_id = eNum;
2395  _edge->setPartitions(partitions);
2396  // put edge number into queue of shared edges to resolve remote partitions
2397  sharedEdgesQueue.push_back(eNum);
2398  }
2399  }
2400  }
2401  }
2402 
2403 #endif
2404 
2405  _nodes.at(1) = nodes.at(inode);
2406  _nodes.at(2) = nodes.at(jnode);
2407  _nodes.at(3) = nodes.at(knode);
2408  _nodes.at(4) = irregular_nodes.at(iiedge);
2409  childNum = mesh->giveNumberOfElements() + 1;
2410  child = new Subdivision :: RS_Tetra(childNum, mesh, this->number, _nodes); // number, parent, coords
2411  mesh->addElement(child);
2412  children.at(1) = childNum;
2413  // set neigbour info
2414  if ( irregulars1 ) {
2415  child->setNeighbor( 1, this->giveNeighbor(nside) );
2416  child->setNeighbor( 2, this->giveNeighbor(kside) );
2417  child->setNeighbor( 4, this->giveNeighbor(jside) );
2418 
2419  child->setIrregular( 1, irregular_nodes.at(kedge) );
2420  child->setIrregular( 3, irregular_nodes.at(jedge) );
2421  child->setIrregular( 2, irregular_nodes.at(iedge) );
2422  } else {
2423  child->setNeighbor( 1, -this->giveNeighbor(nside) );
2424  child->setNeighbor( 2, -this->giveNeighbor(kside) );
2425  child->setNeighbor( 4, -this->giveNeighbor(jside) );
2426  }
2427 
2428  // neihgbor3 of child1 is changed to negative during subdivision (if any) of child2
2429  child->setNeighbor(3, childNum + 1);
2430 #ifdef __PARALLEL_MODE
2431  // set shared info
2432  if ( ishared || jshared || kshared || iishared ) {
2433  child->makeSharedEdges();
2434  if ( ishared ) {
2435  child->setSharedEdge( 2, shared_edges.at(iedge) );
2436  }
2437 
2438  if ( jshared ) {
2439  child->setSharedEdge( 3, shared_edges.at(jedge) );
2440  }
2441 
2442  if ( kshared ) {
2443  child->setSharedEdge( 1, shared_edges.at(kedge) );
2444  }
2445 
2446  if ( iishared ) {
2447  child->setSharedEdge( 4, shared_edges.at(iiedge) );
2448  }
2449 
2450  child->setSharedEdge(5, j_shared_id);
2451  child->setSharedEdge(6, k_shared_id);
2452  }
2453 
2454 #endif
2455 
2456 #ifdef DEBUG_INFO
2457  #ifdef __PARALLEL_MODE
2458  OOFEM_LOG_INFO("[%d] Child %d generated on parent %d [%d] (leIndex %d, nds %d %d %d %d [%d %d %d %d], ngbs %d %d %d %d, irr %d %d %d %d %d %d [%d %d %d %d %d %d])\n",
2459  mesh->giveSubdivision()->giveRank(), childNum,
2460  this->number, this->giveGlobalNumber(), this->leIndex,
2461  _nodes.at(1), _nodes.at(2), _nodes.at(3), _nodes.at(4),
2462  mesh->giveNode( _nodes.at(1) )->giveGlobalNumber(),
2463  mesh->giveNode( _nodes.at(2) )->giveGlobalNumber(),
2464  mesh->giveNode( _nodes.at(3) )->giveGlobalNumber(),
2465  mesh->giveNode( _nodes.at(4) )->giveGlobalNumber(),
2466  child->giveNeighbor(1), child->giveNeighbor(2), child->giveNeighbor(3), child->giveNeighbor(4),
2467  child->giveIrregular(1), child->giveIrregular(2), child->giveIrregular(3),
2468  child->giveIrregular(4), child->giveIrregular(5), child->giveIrregular(6),
2469  ( child->giveIrregular(1) ) ? mesh->giveNode( child->giveIrregular(1) )->giveGlobalNumber() : 0,
2470  ( child->giveIrregular(2) ) ? mesh->giveNode( child->giveIrregular(2) )->giveGlobalNumber() : 0,
2471  ( child->giveIrregular(3) ) ? mesh->giveNode( child->giveIrregular(3) )->giveGlobalNumber() : 0,
2472  ( child->giveIrregular(4) ) ? mesh->giveNode( child->giveIrregular(4) )->giveGlobalNumber() : 0,
2473  ( child->giveIrregular(5) ) ? mesh->giveNode( child->giveIrregular(5) )->giveGlobalNumber() : 0,
2474  ( child->giveIrregular(6) ) ? mesh->giveNode( child->giveIrregular(6) )->giveGlobalNumber() : 0);
2475  #else
2476  OOFEM_LOG_INFO( "Child %d generated on parent %d (leIndex %d, nds %d %d %d %d, ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
2477  childNum, this->number, this->leIndex, _nodes.at(1), _nodes.at(2), _nodes.at(3), _nodes.at(4),
2478  child->giveNeighbor(1), child->giveNeighbor(2), child->giveNeighbor(3), child->giveNeighbor(4),
2479  child->giveIrregular(1), child->giveIrregular(2), child->giveIrregular(3),
2480  child->giveIrregular(4), child->giveIrregular(5), child->giveIrregular(6) );
2481  #endif
2482 #endif
2483 
2484  _nodes.at(1) = irregular_nodes.at(iiedge);
2485  _nodes.at(2) = nodes.at(jnode);
2486  _nodes.at(3) = nodes.at(knode);
2487  _nodes.at(4) = nodes.at(nnode);
2488  childNum = mesh->giveNumberOfElements() + 1;
2489  child = new Subdivision :: RS_Tetra(childNum, mesh, this->number, _nodes); // number, parent, coords
2490  mesh->addElement(child);
2491  children.at(2) = childNum;
2492 
2493  // set neigbour info
2494  if ( irregulars2 ) {
2495  child->setNeighbor( 2, this->giveNeighbor(kside) );
2496  child->setNeighbor( 3, this->giveNeighbor(iside) );
2497  child->setNeighbor( 4, this->giveNeighbor(jside) );
2498 
2499  child->setIrregular( 5, irregular_nodes.at(jjedge) );
2500  child->setIrregular( 6, irregular_nodes.at(kkedge) );
2501  child->setIrregular( 2, irregular_nodes.at(iedge) );
2502  } else {
2503  child->setNeighbor( 2, -this->giveNeighbor(kside) );
2504  child->setNeighbor( 3, -this->giveNeighbor(iside) );
2505  child->setNeighbor( 4, -this->giveNeighbor(jside) );
2506  }
2507 
2508  // neihgbor1 of child2 is changed to negative during subdivision (if any) of child1
2509  child->setNeighbor(1, childNum - 1);
2510 #ifdef __PARALLEL_MODE
2511  // set shared info
2512  if ( ishared || jjshared || kkshared || iishared ) {
2513  child->makeSharedEdges();
2514  if ( ishared ) {
2515  child->setSharedEdge( 2, shared_edges.at(iedge) );
2516  }
2517 
2518  if ( jjshared ) {
2519  child->setSharedEdge( 5, shared_edges.at(jjedge) );
2520  }
2521 
2522  if ( kkshared ) {
2523  child->setSharedEdge( 6, shared_edges.at(kkedge) );
2524  }
2525 
2526  if ( iishared ) {
2527  child->setSharedEdge( 4, shared_edges.at(iiedge) );
2528  }
2529 
2530  child->setSharedEdge(1, j_shared_id);
2531  child->setSharedEdge(3, k_shared_id);
2532  }
2533 
2534 #endif
2535 
2536 #ifdef __PARALLEL_MODE
2537  // set connected elements
2538  if ( mesh->giveNode( nodes.at(jnode) )->giveParallelMode() == DofManager_shared ) {
2539  mesh->giveNode( nodes.at(jnode) )->eraseConnectedElement( this->giveNumber() );
2540  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum);
2541  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum - 1);
2542  }
2543 
2544  if ( mesh->giveNode( nodes.at(knode) )->giveParallelMode() == DofManager_shared ) {
2545  mesh->giveNode( nodes.at(knode) )->eraseConnectedElement( this->giveNumber() );
2546  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum);
2547  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum - 1);
2548  }
2549 
2550  if ( mesh->giveNode( nodes.at(inode) )->giveParallelMode() == DofManager_shared ) {
2551  mesh->giveNode( nodes.at(inode) )->eraseConnectedElement( this->giveNumber() );
2552  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 1);
2553  }
2554 
2555  if ( mesh->giveNode( nodes.at(nnode) )->giveParallelMode() == DofManager_shared ) {
2556  mesh->giveNode( nodes.at(nnode) )->eraseConnectedElement( this->giveNumber() );
2557  mesh->giveNode( nodes.at(nnode) )->insertConnectedElement(childNum);
2558  }
2559 
2560  if ( iishared ) {
2561  mesh->giveNode( irregular_nodes.at(iiedge) )->insertConnectedElement(childNum);
2562  mesh->giveNode( irregular_nodes.at(iiedge) )->insertConnectedElement(childNum - 1);
2563  }
2564 
2565 #endif
2566 
2567 #ifdef DEBUG_INFO
2568  #ifdef __PARALLEL_MODE
2569  OOFEM_LOG_INFO("[%d] Child %d generated on parent %d [%d] (leIndex %d, nds %d %d %d %d [%d %d %d %d], ngbs %d %d %d %d, irr %d %d %d %d %d %d [%d %d %d %d %d %d])\n",
2570  mesh->giveSubdivision()->giveRank(), childNum,
2571  this->number, this->giveGlobalNumber(), this->leIndex,
2572  _nodes.at(1), _nodes.at(2), _nodes.at(3), _nodes.at(4),
2573  mesh->giveNode( _nodes.at(1) )->giveGlobalNumber(),
2574  mesh->giveNode( _nodes.at(2) )->giveGlobalNumber(),
2575  mesh->giveNode( _nodes.at(3) )->giveGlobalNumber(),
2576  mesh->giveNode( _nodes.at(4) )->giveGlobalNumber(),
2577  child->giveNeighbor(1), child->giveNeighbor(2), child->giveNeighbor(3), child->giveNeighbor(4),
2578  child->giveIrregular(1), child->giveIrregular(2), child->giveIrregular(3),
2579  child->giveIrregular(4), child->giveIrregular(5), child->giveIrregular(6),
2580  ( child->giveIrregular(1) ) ? mesh->giveNode( child->giveIrregular(1) )->giveGlobalNumber() : 0,
2581  ( child->giveIrregular(2) ) ? mesh->giveNode( child->giveIrregular(2) )->giveGlobalNumber() : 0,
2582  ( child->giveIrregular(3) ) ? mesh->giveNode( child->giveIrregular(3) )->giveGlobalNumber() : 0,
2583  ( child->giveIrregular(4) ) ? mesh->giveNode( child->giveIrregular(4) )->giveGlobalNumber() : 0,
2584  ( child->giveIrregular(5) ) ? mesh->giveNode( child->giveIrregular(5) )->giveGlobalNumber() : 0,
2585  ( child->giveIrregular(6) ) ? mesh->giveNode( child->giveIrregular(6) )->giveGlobalNumber() : 0);
2586  #else
2587  OOFEM_LOG_INFO( "Child %d generated on parent %d (leIndex %d, nds %d %d %d %d, ngbs %d %d %d %d, irr %d %d %d %d %d %d)\n",
2588  childNum, this->number, this->leIndex, _nodes.at(1), _nodes.at(2), _nodes.at(3), _nodes.at(4),
2589  child->giveNeighbor(1), child->giveNeighbor(2), child->giveNeighbor(3), child->giveNeighbor(4),
2590  child->giveIrregular(1), child->giveIrregular(2), child->giveIrregular(3),
2591  child->giveIrregular(4), child->giveIrregular(5), child->giveIrregular(6) );
2592  #endif
2593 #endif
2594 
2595  // set connected elements to nonshared nodes on outer boundary
2596 #ifdef __PARALLEL_MODE
2597  if ( !iishared ) {
2598 #endif
2599  if ( mesh->giveNode( irregular_nodes.at(iiedge) )->giveNumber() < 0 ) { // check for marked local irregular
2600  // irregular node on outer boundary
2601  // insert connected elements
2602  mesh->giveNode( irregular_nodes.at(iiedge) )->insertConnectedElement(childNum);
2603  mesh->giveNode( irregular_nodes.at(iiedge) )->insertConnectedElement(childNum - 1);
2604  }
2605 
2606  // update connectivity of both end nodes of iiedge if not shared
2607  // and ONLY if there is already list of connected elements
2608 #ifdef __PARALLEL_MODE
2609  if ( mesh->giveNode( nodes.at(inode) )->giveParallelMode() != DofManager_shared ) { // update already done
2610 #endif
2611  if ( mesh->giveNode( nodes.at(inode) )->giveConnectedElements()->giveSize() ) {
2612  mesh->giveNode( nodes.at(inode) )->eraseConnectedElement( this->giveNumber() );
2613  mesh->giveNode( nodes.at(inode) )->insertConnectedElement(childNum - 1);
2614  }
2615 
2616 #ifdef __PARALLEL_MODE
2617  }
2618 #endif
2619 #ifdef __PARALLEL_MODE
2620  if ( mesh->giveNode( nodes.at(nnode) )->giveParallelMode() != DofManager_shared ) { // update already done
2621 #endif
2622  if ( mesh->giveNode( nodes.at(nnode) )->giveConnectedElements()->giveSize() ) {
2623  mesh->giveNode( nodes.at(nnode) )->eraseConnectedElement( this->giveNumber() );
2624  mesh->giveNode( nodes.at(nnode) )->insertConnectedElement(childNum);
2625  }
2626 
2627 #ifdef __PARALLEL_MODE
2628  }
2629 #endif
2630 #ifdef __PARALLEL_MODE
2631  }
2632 #endif
2633 #ifdef __PARALLEL_MODE
2634  if ( !ishared ) {
2635 #endif
2636  // update connectivity of both end nodes of iedge if not shared
2637  // and ONLY if there is already list of connected elements
2638 #ifdef __PARALLEL_MODE
2639  if ( mesh->giveNode( nodes.at(jnode) )->giveParallelMode() != DofManager_shared ) { // update already done
2640 #endif
2641  if ( mesh->giveNode( nodes.at(jnode) )->giveConnectedElements()->giveSize() ) {
2642  mesh->giveNode( nodes.at(jnode) )->eraseConnectedElement( this->giveNumber() );
2643  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum);
2644  mesh->giveNode( nodes.at(jnode) )->insertConnectedElement(childNum - 1);
2645  }
2646 
2647 #ifdef __PARALLEL_MODE
2648  }
2649 #endif
2650 #ifdef __PARALLEL_MODE
2651  if ( mesh->giveNode( nodes.at(knode) )->giveParallelMode() != DofManager_shared ) { // update already done
2652 #endif
2653  if ( mesh->giveNode( nodes.at(knode) )->giveConnectedElements()->giveSize() ) {
2654  mesh->giveNode( nodes.at(knode) )->eraseConnectedElement( this->giveNumber() );
2655  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum);
2656  mesh->giveNode( nodes.at(knode) )->insertConnectedElement(childNum - 1);
2657  }
2658 
2659 #ifdef __PARALLEL_MODE
2660  }
2661 #endif
2662 #ifdef __PARALLEL_MODE
2663  }
2664 #endif
2665  }
2666 
2667  if ( irregulars1 ) {
2668 #ifdef DEBUG_CHECK
2669  if ( !mesh->giveElement( children.at(1) )->giveIrregular(leIndex1) ) {
2670  OOFEM_ERROR("leIndex incosistency for child 1 of %d", this->number);
2671  }
2672 
2673 #endif
2674  mesh->giveElement( children.at(1) )->setLeIndex(leIndex1);
2675  mesh->giveElement( children.at(1) )->generate(sharedEdgesQueue);
2676  }
2677 
2678  if ( irregulars2 ) {
2679 #ifdef DEBUG_CHECK
2680  if ( !mesh->giveElement( children.at(2) )->giveIrregular(leIndex2) ) {
2681  OOFEM_ERROR("leIndex incosistency for child 2 of %d", this->number);
2682  }
2683 
2684 #endif
2685  mesh->giveElement( children.at(2) )->setLeIndex(leIndex2);
2686  mesh->giveElement( children.at(2) )->generate(sharedEdgesQueue);
2687  }
2688 
2689  // if there is neighbor of "this" (local or remote) not designated for bisection
2690  // its neihgbor corresponding to "this" must be made negative to enforce update_neighbors
2691  for ( i = 1; i <= 4; i++ ) {
2692  if ( this->neghbours_base_elements.at(i) ) {
2693 #ifdef DEBUG_CHECK
2694  if ( this->neghbours_base_elements.at(i) < 0 ) {
2695  OOFEM_ERROR("negative neighbor of %d not expected", this->number);
2696  }
2697 
2698 #endif
2699 
2700  ngb = mesh->giveElement( this->neghbours_base_elements.at(i) );
2701  if ( !ngb->hasIrregulars() ) {
2702  ind = ngb->giveNeighbors()->findFirstIndexOf(this->number);
2703  if ( ind ) {
2704  ngb->setNeighbor(ind, -this->number);
2705  }
2706  }
2707  }
2708  }
2709  }
2710 }
2711 
2712 
2713 
2714 void
2716 {
2717  bool found;
2718  int i, iside, parentNeighbor;
2719  const IntArray *ca;
2720 
2721  /*
2722  * neghbours_base_elements pre-set in generate using this rule:
2723  * value > 0 .. real neighbour known, value equals to real neighbor on new mesh
2724  * value = 0 .. boundary
2725  * value < 0 .. neighbour set to parent element, actual neighbour has to be determined from parent children
2726  * or if no child exists then to parent new counterpart.
2727  */
2728 
2729  for ( iside = 1; iside <= 3; iside++ ) {
2730  if ( this->neghbours_base_elements.at(iside) < 0 ) {
2731  parentNeighbor = -this->neghbours_base_elements.at(iside);
2732  // test if parentNeighbor has been subdivided
2733  if ( mesh->giveElement(parentNeighbor)->hasIrregulars() ) { // HUHU replace by !isTerminal
2734  // old neighbour bisected -> we loop over its children to find an appropriate new neighbor
2735  ca = mesh->giveElement(parentNeighbor)->giveChildren();
2736  found = false;
2737  for ( i = 1; i <= ca->giveSize(); i++ ) {
2738  if ( this->isNeighborOf( mesh->giveElement( ca->at(i) ) ) ) {
2739  this->neghbours_base_elements.at(iside) = ca->at(i);
2740  found = true;
2741  break;
2742  }
2743  }
2744 
2745  if ( !found ) {
2746  OOFEM_ERROR("update_neighbours failed for element %d", this->number);
2747  }
2748  } else {
2749  // parent neighbour remains actual neighbour
2750  this->neghbours_base_elements.at(iside) = parentNeighbor;
2751  }
2752  } else if ( this->neghbours_base_elements.at(iside) > 0 ) {
2753  // neighbor element already set
2754  }
2755  } // end loop over element sides
2756 }
2757 
2758 
2759 void
2761 {
2762  bool found;
2763  int i, j, k, iside, parentNeighbor;
2764  const IntArray *ca1, *ca2, *ca3;
2765 
2766  /*
2767  * neghbours_base_elements pre-set in generate using this rule:
2768  * value > 0 .. real neighbour known, value equals to real neighbor on new mesh
2769  * value = 0 .. boundary
2770  * value < 0 .. neighbour set to parent element, actual neighbour has to be determined from terminal children of parent
2771  * or if no child exists then to parent new counterpart.
2772  */
2773 
2774  for ( iside = 1; iside <= 4; iside++ ) {
2775  if ( this->neghbours_base_elements.at(iside) < 0 ) {
2776  parentNeighbor = -this->neghbours_base_elements.at(iside);
2777  // test if parentNeighbor has been subdivided
2778  if ( mesh->giveElement(parentNeighbor)->hasIrregulars() ) { // HUHU replace by !isTerminal
2779  // old neighbour bisected -> we loop over its terminal children to find an appropriate new neighbor;
2780  // there may be at maximum 3 levels of parent tetra subdivision
2781  ca1 = mesh->giveElement(parentNeighbor)->giveChildren();
2782  found = false;
2783  for ( i = 1; i <= ca1->giveSize(); i++ ) {
2784  if ( mesh->giveElement( ca1->at(i) )->isTerminal() ) {
2785  if ( this->isNeighborOf( mesh->giveElement( ca1->at(i) ) ) ) {
2786  this->neghbours_base_elements.at(iside) = ca1->at(i);
2787  found = true;
2788  break;
2789  }
2790  } else {
2791  ca2 = mesh->giveElement( ca1->at(i) )->giveChildren();
2792  for ( j = 1; j <= ca2->giveSize(); j++ ) {
2793  if ( mesh->giveElement( ca2->at(j) )->isTerminal() ) {
2794  if ( this->isNeighborOf( mesh->giveElement( ca2->at(j) ) ) ) {
2795  this->neghbours_base_elements.at(iside) = ca2->at(j);
2796  found = true;
2797  break;
2798  }
2799  } else {
2800  ca3 = mesh->giveElement( ca2->at(j) )->giveChildren();
2801  for ( k = 1; k <= ca3->giveSize(); k++ ) {
2802  if ( mesh->giveElement( ca3->at(k) )->isTerminal() ) {
2803  if ( this->isNeighborOf( mesh->giveElement( ca3->at(k) ) ) ) {
2804  this->neghbours_base_elements.at(iside) = ca3->at(k);
2805  found = true;
2806  break;
2807  }
2808  }
2809 
2810 #ifdef DEBUG_CHECK
2811  else {
2812  OOFEM_ERROR("too many subdivision levels for element %d", this->number);
2813  }
2814 #endif
2815  }
2816  }
2817 
2818  if ( found ) {
2819  break;
2820  }
2821  }
2822  }
2823 
2824  if ( found ) {
2825  break;
2826  }
2827  }
2828 
2829  if ( !found ) {
2830  OOFEM_ERROR("failed for element %d (side %d, elem %d)",
2831  this->number, iside, parentNeighbor);
2832  }
2833  } else {
2834  // parent neighbour remains actual neighbour
2835  this->neghbours_base_elements.at(iside) = parentNeighbor;
2836  }
2837  } else if ( this->neghbours_base_elements.at(iside) > 0 ) {
2838  // neighbor element already set
2839  }
2840  } // end loop over element side faces
2841 
2842 #ifdef DEBUG_CHECK
2843  // check updated neighbors
2844  IntArray snodes1, snodes2;
2845  RS_Element *ngb;
2846  for ( i = 1; i <= 4; i++ ) {
2847  if ( this->neghbours_base_elements.at(i) ) {
2848  if ( this->neghbours_base_elements.at(i) > 0 ) {
2849  this->giveSideNodes(i, snodes1);
2850  ngb = mesh->giveElement( this->neghbours_base_elements.at(i) );
2851  if ( ngb->giveNumber() < 0 ) {
2852  OOFEM_ERROR("neighbour %d of %d is temporary",
2853  this->neghbours_base_elements.at(i), this->number);
2854  }
2855 
2856  if ( !ngb->isTerminal() ) {
2857  OOFEM_ERROR("neighbour %d of %d not terminal",
2858  this->neghbours_base_elements.at(i), this->number);
2859  }
2860 
2861  if ( ngb->giveNumber() < this->number ) {
2862  // ngb has been already processed
2863  j = ngb->giveNeighbors()->findFirstIndexOf(this->number);
2864  if ( j ) {
2865  ngb->giveSideNodes(j, snodes2);
2866  for ( k = 1; k <= 3; k++ ) {
2867  if ( snodes2.findFirstIndexOf( snodes1.at(k) ) ) {
2868  continue;
2869  }
2870 
2871  OOFEM_ERROR( "element %d is not neighbor (%d) of element %d",
2872  this->number, i, this->neghbours_base_elements.at(i) );
2873  }
2874  } else {
2875  OOFEM_ERROR( "element %d is not neighbor (%d) of element %d",
2876  this->number, i, this->neghbours_base_elements.at(i) );
2877  }
2878  } else {
2879  // ngb has not been yet processed ==> I cannot check particular side of ngb
2880  for ( k = 1; k <= 3; k++ ) {
2881  if ( ngb->giveNodes()->findFirstIndexOf( snodes1.at(k) ) ) {
2882  continue;
2883  }
2884 
2885  OOFEM_ERROR( "element %d is not neighbor of element %d",
2886  this->number, this->neghbours_base_elements.at(i) );
2887  }
2888  }
2889  } else {
2890  OOFEM_ERROR("negative neighbour %d of %d not expected",
2891  this->neghbours_base_elements.at(i), this->number);
2892  }
2893  }
2894  }
2895 
2896 #endif
2897 }
2898 
2899 
2900 bool
2902 {
2903  // Simplified implementation, considering only one element type - triangles
2904  int i, _c = 0;
2905  for ( i = 1; i <= 3; i++ ) {
2906  _c += elem->containsNode( this->nodes.at(i) );
2907  }
2908 
2909  return ( _c == 2 );
2910 }
2911 
2912 
2913 bool
2915 {
2916  // Simplified implementation, considering only one element type - tetras
2917  int i, _c = 0;
2918  for ( i = 1; i <= 4; i++ ) {
2919  _c += elem->containsNode( this->nodes.at(i) );
2920  }
2921 
2922  return ( _c == 3 );
2923 }
2924 
2925 
2926 void
2928 {
2929  int inode, jnode;
2930 
2931  inode = iside;
2932  jnode = ( iside < 3 ) ? iside + 1 : 1;
2933 
2934  snodes.resize(2);
2935  snodes.at(1) = nodes.at(inode);
2936  snodes.at(2) = nodes.at(jnode);
2937 }
2938 
2939 
2940 void
2942 {
2943  int inode, jnode, knode;
2944 
2945  if ( iside == 1 ) {
2946  inode = 1;
2947  jnode = 2;
2948  knode = 3;
2949  } else {
2950  inode = iside - 1;
2951  jnode = ( iside < 4 ) ? iside : 1;
2952  knode = 4;
2953  }
2954 
2955  snodes.resize(3);
2956  snodes.at(1) = nodes.at(inode);
2957  snodes.at(2) = nodes.at(jnode);
2958  snodes.at(3) = nodes.at(knode);
2959 }
2960 
2961 
2962 
2963 double
2965 {
2966  double x1, x2, x3, y1, y2, y3;
2967  x1 = mesh->giveNode( nodes.at(1) )->giveCoordinate(1);
2968  x2 = mesh->giveNode( nodes.at(2) )->giveCoordinate(1);
2969  x3 = mesh->giveNode( nodes.at(3) )->giveCoordinate(1);
2970 
2971  y1 = mesh->giveNode( nodes.at(1) )->giveCoordinate(2);
2972  y2 = mesh->giveNode( nodes.at(2) )->giveCoordinate(2);
2973  y3 = mesh->giveNode( nodes.at(3) )->giveCoordinate(2);
2974 
2975  return sqrt( fabs( 0.5 * ( x2 * y3 + x1 * y2 + y1 * x3 - x2 * y1 - x3 * y2 - x1 * y3 ) ) );
2976 }
2977 
2978 
2979 double
2981 {
2982  double x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4;
2983  double dx1, dx2, dx3, dy1, dy2, dy3, dz1, dz2, dz3;
2984  double vol, d;
2985  x1 = mesh->giveNode( nodes.at(1) )->giveCoordinate(1);
2986  x2 = mesh->giveNode( nodes.at(2) )->giveCoordinate(1);
2987  x3 = mesh->giveNode( nodes.at(3) )->giveCoordinate(1);
2988  x4 = mesh->giveNode( nodes.at(4) )->giveCoordinate(1);
2989 
2990  y1 = mesh->giveNode( nodes.at(1) )->giveCoordinate(2);
2991  y2 = mesh->giveNode( nodes.at(2) )->giveCoordinate(2);
2992  y3 = mesh->giveNode( nodes.at(3) )->giveCoordinate(2);
2993  y4 = mesh->giveNode( nodes.at(4) )->giveCoordinate(2);
2994 
2995  z1 = mesh->giveNode( nodes.at(1) )->giveCoordinate(3);
2996  z2 = mesh->giveNode( nodes.at(2) )->giveCoordinate(3);
2997  z3 = mesh->giveNode( nodes.at(3) )->giveCoordinate(3);
2998  z4 = mesh->giveNode( nodes.at(4) )->giveCoordinate(3);
2999 
3000  dx1 = x2 - x1;
3001  dy1 = y2 - y1;
3002  dz1 = z2 - z1;
3003  dx2 = x3 - x1;
3004  dy2 = y3 - y1;
3005  dz2 = z3 - z1;
3006  dx3 = x4 - x1;
3007  dy3 = y4 - y1;
3008  dz3 = z4 - z1;
3009 
3010  vol = ( dx3 * ( dy1 * dz2 - dz1 * dy2 ) + dy3 * ( dz1 * dx2 - dx1 * dz2 ) + dz3 * ( dx1 * dy2 - dy1 * dx2 ) ) / 6.0;
3011  d = exp(log(vol) / 3.0);
3012 
3013  return d;
3014 }
3015 
3016 
3017 void
3019 {
3020  IntArray me(1), conn;
3021  int iNode, jNode, el, i;
3022 
3024  neghbours_base_elements.zero(); // initialized to have no neighbour
3025 
3026 #if 0 // old version
3027  me.at(1) = this->number;
3028  ct->giveElementNeighbourList(conn, me);
3029 
3030  int iside;
3031  for ( iside = 1; iside <= 3; iside++ ) {
3032  iNode = nodes.at(iside);
3033  jNode = nodes.at( ( iside == 3 ) ? 1 : iside + 1 );
3034 
3035  // select right neighbour
3036  for ( i = 1; i <= conn.giveSize(); i++ ) {
3037  el = conn.at(i);
3038  if ( el == this->number ) {
3039  continue;
3040  }
3041 
3042  #ifdef __PARALLEL_MODE
3043  if ( mesh->giveElement(el)->giveParallelMode() != Element_local ) {
3044  continue;
3045  }
3046 
3047  #endif
3048  if ( mesh->giveElement(el)->containsNode(iNode) &&
3049  mesh->giveElement(el)->containsNode(jNode) ) {
3050  neghbours_base_elements.at(iside) = el;
3051  break;
3052  }
3053  }
3054  }
3055 
3056 #endif
3057 
3058  me.at(1) = nodes.at(3);
3059  ct->giveNodeNeighbourList(conn, me);
3060  iNode = nodes.at(1);
3061  jNode = nodes.at(2);
3062 
3063  // select right neighbour
3064  for ( i = 1; i <= conn.giveSize(); i++ ) {
3065  el = conn.at(i);
3066  if ( el == this->number ) {
3067  continue;
3068  }
3069 
3070 #ifdef __PARALLEL_MODE
3071  if ( mesh->giveElement(el)->giveParallelMode() != Element_local ) {
3072  continue;
3073  }
3074 
3075 #endif
3076 
3077  if ( mesh->giveElement(el)->containsNode(iNode) ) {
3078  neghbours_base_elements.at(3) = el;
3079  break;
3080  }
3081  }
3082 
3083  for ( i = 1; i <= conn.giveSize(); i++ ) {
3084  el = conn.at(i);
3085  if ( el == this->number ) {
3086  continue;
3087  }
3088 
3089 #ifdef __PARALLEL_MODE
3090  if ( mesh->giveElement(el)->giveParallelMode() != Element_local ) {
3091  continue;
3092  }
3093 
3094 #endif
3095 
3096  if ( mesh->giveElement(el)->containsNode(jNode) ) {
3097  neghbours_base_elements.at(2) = el;
3098  break;
3099  }
3100  }
3101 
3102  me.at(1) = iNode;
3103  ct->giveNodeNeighbourList(conn, me);
3104 
3105  // select right neighbour
3106  for ( i = 1; i <= conn.giveSize(); i++ ) {
3107  el = conn.at(i);
3108  if ( el == this->number ) {
3109  continue;
3110  }
3111 
3112 #ifdef __PARALLEL_MODE
3113  if ( mesh->giveElement(el)->giveParallelMode() != Element_local ) {
3114  continue;
3115  }
3116 
3117 #endif
3118 
3119  if ( mesh->giveElement(el)->containsNode(jNode) ) {
3120  neghbours_base_elements.at(1) = el;
3121  break;
3122  }
3123  }
3124 }
3125 
3126 
3127 void
3129 {
3130  IntArray me(1), conn;
3131  int iNode, jNode, kNode, el, i;
3132 
3134  neghbours_base_elements.zero(); // initialized to have no neighbour
3135 
3136 #if 0 // old version
3137  me.at(1) = this->number;
3138  ct->giveElementNeighbourList(conn, me);
3139 
3140  iNode = nodes.at(1);
3141  jNode = nodes.at(2);
3142  kNode = nodes.at(3);
3143 
3144  // select right base neighbour
3145  for ( i = 1; i <= conn.giveSize(); i++ ) {
3146  el = conn.at(i);
3147  if ( el == this->number ) {
3148  continue;
3149  }
3150 
3151  #ifdef __PARALLEL_MODE
3152  if ( mesh->giveElement(el)->giveParallelMode() != Element_local ) {
3153  continue;
3154  }
3155 
3156  #endif
3157  if ( mesh->giveElement(el)->containsNode(iNode) &&
3158  mesh->giveElement(el)->containsNode(jNode) &&
3159  mesh->giveElement(el)->containsNode(kNode) ) {
3160  neghbours_base_elements.at(1) = el;
3161  break;
3162  }
3163  }
3164 
3165  kNode = nodes.at(4);
3166  int iside;
3167  for ( iside = 1; iside <= 3; iside++ ) {
3168  iNode = nodes.at(iside);
3169  jNode = nodes.at( ( iside == 3 ) ? 1 : iside + 1 );
3170 
3171  // select right side neighbour
3172  for ( i = 1; i <= conn.giveSize(); i++ ) {
3173  el = conn.at(i);
3174  if ( el == this->number ) {
3175  continue;
3176  }
3177 
3178  #ifdef __PARALLEL_MODE
3179  if ( mesh->giveElement(el)->giveParallelMode() != Element_local ) {
3180  continue;
3181  }
3182 
3183  #endif
3184  if ( mesh->giveElement(el)->containsNode(iNode) &&
3185  mesh->giveElement(el)->containsNode(jNode) &&
3186  mesh->giveElement(el)->containsNode(kNode) ) {
3187  neghbours_base_elements.at(iside + 1) = el;
3188  break;
3189  }
3190  }
3191  }
3192 
3193 #endif
3194 
3195  bool has_i, has_j, has_k;
3196 
3197  me.at(1) = nodes.at(4);
3198  ct->giveNodeNeighbourList(conn, me);
3199  iNode = nodes.at(1);
3200  jNode = nodes.at(2);
3201  kNode = nodes.at(3);
3202 
3203  // select right neighbour
3204  for ( i = 1; i <= conn.giveSize(); i++ ) {
3205  el = conn.at(i);
3206  if ( el == this->number ) {
3207  continue;
3208  }
3209 
3210 #ifdef __PARALLEL_MODE
3211  if ( mesh->giveElement(el)->giveParallelMode() != Element_local ) {
3212  continue;
3213  }
3214 
3215 #endif
3216 
3217  has_i = mesh->giveElement(el)->containsNode(iNode);
3218  has_j = mesh->giveElement(el)->containsNode(jNode);
3219  has_k = mesh->giveElement(el)->containsNode(kNode);
3220 
3221  if ( has_i && has_j ) {
3222  neghbours_base_elements.at(2) = el;
3223  }
3224 
3225  if ( has_j && has_k ) {
3226  neghbours_base_elements.at(3) = el;
3227  }
3228 
3229  if ( has_k && has_i ) {
3230  neghbours_base_elements.at(4) = el;
3231  }
3232  }
3233 
3234  me.at(1) = iNode;
3235  ct->giveNodeNeighbourList(conn, me);
3236 
3237  // select right neighbour
3238  for ( i = 1; i <= conn.giveSize(); i++ ) {
3239  el = conn.at(i);
3240  if ( el == this->number ) {
3241  continue;
3242  }
3243 
3244 #ifdef __PARALLEL_MODE
3245  if ( mesh->giveElement(el)->giveParallelMode() != Element_local ) {
3246  continue;
3247  }
3248 
3249 #endif
3250 
3251  if ( mesh->giveElement(el)->containsNode(jNode) &&
3252  mesh->giveElement(el)->containsNode(kNode) ) {
3253  neghbours_base_elements.at(1) = el;
3254  break;
3255  }
3256  }
3257 
3258  /*
3259  * #ifdef DEBUG_INFO
3260  * OOFEM_LOG_INFO("Elem %d connectivity imported (nds %d %d %d %d ngbs %d %d %d %d)\n", this->number,
3261  * nodes.at(1), nodes.at(2), nodes.at(3), nodes.at(4),
3262  * neghbours_base_elements.at(1), neghbours_base_elements.at(2),
3263  * neghbours_base_elements.at(3), neghbours_base_elements.at(4));
3264  ******#endif
3265  */
3266 }
3267 
3268 
3269 #ifdef __OOFEG
3270 void
3272 //
3273 // draws graphics representation of receiver
3274 //
3275 {
3276  GraphicObj *go;
3277  EPixel color;
3278  BOOLEAN suc;
3279  const char *colors[] = {
3280  "orange", "black"
3281  };
3282 
3283  WCRec p [ 1 ]; /* point */
3284  p [ 0 ].x = ( FPNum ) this->giveCoordinate(1);
3285  p [ 0 ].y = ( FPNum ) this->giveCoordinate(2);
3286  p [ 0 ].z = ( FPNum ) this->giveCoordinate(3);
3287 
3288  EASValsSetMType(FILLED_CIRCLE_MARKER);
3289  color = ColorGetPixelFromString(const_cast< char * >(colors [ 0 ]), & suc);
3290  EASValsSetColor(color);
3291  //EASValsSetColor( gc.getNodeColor() );
3292  EASValsSetLayer(OOFEG_RAW_GEOMETRY_LAYER);
3293  EASValsSetMSize(8);
3294  go = CreateMarker3D(p);
3295  EGWithMaskChangeAttributes(COLOR_MASK | LAYER_MASK | MTYPE_MASK | MSIZE_MASK, go);
3296  EMAddGraphicsToModel(ESIModel(), go);
3297 
3298  char num [ 6 ];
3299  color = ColorGetPixelFromString(const_cast< char * >(colors [ 1 ]), & suc);
3300  EASValsSetColor(color);
3301  //EASValsSetColor( gc.getNodeColor() );
3302  EASValsSetLayer(OOFEG_RAW_GEOMETRY_LAYER);
3303  p [ 0 ].x = ( FPNum ) this->giveCoordinate(1);
3304  p [ 0 ].y = ( FPNum ) this->giveCoordinate(2);
3305  p [ 0 ].z = ( FPNum ) this->giveCoordinate(3);
3306  sprintf(num, "%d", this->number);
3307  go = CreateAnnText3D(p, num);
3308  EGWithMaskChangeAttributes(COLOR_MASK | LAYER_MASK, go);
3309  EMAddGraphicsToModel(ESIModel(), go);
3310 }
3311 
3312 
3313 void
3315 {
3316  WCRec p [ 3 ];
3317  GraphicObj *go;
3318  EPixel color;
3319  BOOLEAN suc;
3320  const char *colors[] = {
3321  "DodgerBlue", "black"
3322  };
3323 
3324  EASValsSetLineWidth(OOFEG_RAW_GEOMETRY_WIDTH);
3325  color = ColorGetPixelFromString(const_cast< char * >(colors [ 0 ]), & suc);
3326  EASValsSetColor(color);
3327  color = ColorGetPixelFromString(const_cast< char * >(colors [ 1 ]), & suc);
3328  EASValsSetEdgeColor(color);
3329  //EASValsSetColor( gc.getElementColor() );
3330  //EASValsSetEdgeColor( gc.getElementEdgeColor() );
3331  EASValsSetEdgeFlag(true);
3332  EASValsSetLayer(OOFEG_RAW_GEOMETRY_LAYER);
3333  // EASValsSetShrink(0.8);
3334  p [ 0 ].x = ( FPNum ) mesh->giveNode( nodes.at(1) )->giveCoordinate(1);
3335  p [ 0 ].y = ( FPNum ) mesh->giveNode( nodes.at(1) )->giveCoordinate(2);
3336  p [ 0 ].z = 0.;
3337  p [ 1 ].x = ( FPNum ) mesh->giveNode( nodes.at(2) )->giveCoordinate(1);
3338  p [ 1 ].y = ( FPNum ) mesh->giveNode( nodes.at(2) )->giveCoordinate(2);
3339  p [ 1 ].z = 0.;
3340  p [ 2 ].x = ( FPNum ) mesh->giveNode( nodes.at(3) )->giveCoordinate(1);
3341  p [ 2 ].y = ( FPNum ) mesh->giveNode( nodes.at(3) )->giveCoordinate(2);
3342  p [ 2 ].z = 0.;
3343 
3344  go = CreateTriangle3D(p);
3345  //EGWithMaskChangeAttributes(WIDTH_MASK | COLOR_MASK | SHRINK_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
3346  EGWithMaskChangeAttributes(WIDTH_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
3347  EGAttachObject(go, ( EObjectP ) this);
3348  EMAddGraphicsToModel(ESIModel(), go);
3349 }
3350 
3351 
3352 void
3354 {
3355  WCRec p [ 4 ];
3356  GraphicObj *go;
3357  EPixel color;
3358  BOOLEAN suc;
3359  const char *colors[] = {
3360  "DodgerBlue", "black"
3361  };
3362 
3363  EASValsSetLineWidth(OOFEG_RAW_GEOMETRY_WIDTH);
3364  color = ColorGetPixelFromString(const_cast< char * >(colors [ 0 ]), & suc);
3365  EASValsSetColor(color);
3366  color = ColorGetPixelFromString(const_cast< char * >(colors [ 1 ]), & suc);
3367  EASValsSetEdgeColor(color);
3368  //EASValsSetColor( gc.getElementColor() );
3369  //EASValsSetEdgeColor( gc.getElementEdgeColor() );
3370  EASValsSetEdgeFlag(true);
3371  EASValsSetLayer(OOFEG_RAW_GEOMETRY_LAYER);
3372  EASValsSetFillStyle(FILL_SOLID);
3373  //EASValsSetShrink(0.8);
3374  p [ 0 ].x = ( FPNum ) mesh->giveNode( nodes.at(1) )->giveCoordinate(1);
3375  p [ 0 ].y = ( FPNum ) mesh->giveNode( nodes.at(1) )->giveCoordinate(2);
3376  p [ 0 ].z = ( FPNum ) mesh->giveNode( nodes.at(1) )->giveCoordinate(3);
3377  p [ 1 ].x = ( FPNum ) mesh->giveNode( nodes.at(2) )->giveCoordinate(1);
3378  p [ 1 ].y = ( FPNum ) mesh->giveNode( nodes.at(2) )->giveCoordinate(2);
3379  p [ 1 ].z = ( FPNum ) mesh->giveNode( nodes.at(2) )->giveCoordinate(3);
3380  p [ 2 ].x = ( FPNum ) mesh->giveNode( nodes.at(3) )->giveCoordinate(1);
3381  p [ 2 ].y = ( FPNum ) mesh->giveNode( nodes.at(3) )->giveCoordinate(2);
3382  p [ 2 ].z = ( FPNum ) mesh->giveNode( nodes.at(3) )->giveCoordinate(3);
3383  p [ 3 ].x = ( FPNum ) mesh->giveNode( nodes.at(4) )->giveCoordinate(1);
3384  p [ 3 ].y = ( FPNum ) mesh->giveNode( nodes.at(4) )->giveCoordinate(2);
3385  p [ 3 ].z = ( FPNum ) mesh->giveNode( nodes.at(4) )->giveCoordinate(3);
3386 
3387  go = CreateTetra(p);
3388  //EGWithMaskChangeAttributes(WIDTH_MASK | COLOR_MASK | SHRINK_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK | FILL_MASK, go);
3389  EGWithMaskChangeAttributes(WIDTH_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK | FILL_MASK, go);
3390  EGAttachObject(go, ( EObjectP ) this);
3391  EMAddGraphicsToModel(ESIModel(), go);
3392 }
3393 #endif
3394 
3395 
3397 Subdivision :: createMesh(TimeStep *tStep, int domainNumber, int domainSerNum, Domain **dNew)
3398 {
3399  // import data from old domain
3400  int parent, nnodes = domain->giveNumberOfDofManagers(), nelems = domain->giveNumberOfElements();
3401  IntArray enodes;
3402  Subdivision :: RS_Node *_node;
3403  Subdivision :: RS_Element *_element;
3404  Timer timer;
3405 
3406  timer.startTimer();
3407 
3408  if ( this->mesh ) {
3409  delete mesh;
3410  }
3411 
3412  mesh = new Subdivision :: RS_Mesh(this);
3413 
3414  // import nodes
3415  // all nodes (including remote which are not needed) are imported to ensure consistency between
3416  // node number (in the mesh) and its parent number (in the domain) because the domain is used to import connectivities
3417  for ( int i = 1; i <= nnodes; i++ ) {
3418  _node = new Subdivision :: RS_Node( i, mesh, i, * ( domain->giveNode ( i )->giveCoordinates() ),
3420  domain->giveNode ( i )->isBoundary() );
3422 #ifdef __PARALLEL_MODE
3424  _node->setPartitions( * domain->giveNode(i)->givePartitionList() );
3425 #endif
3426  this->mesh->addNode(_node);
3427  }
3428 
3429  // import elements
3430  // all elements (including remote which are not needed) are imported to ensure consistency between
3431  // element number (in the mesh) and its parent number (in the domain) because the domain is used to import connectivities
3432  for ( int i = 1; i <= nelems; i++ ) {
3433  if ( domain->giveElement(i)->giveGeometryType() == EGT_triangle_1 ) {
3434  enodes.resize(3);
3435  for ( int j = 1; j <= 3; j++ ) {
3436  enodes.at(j) = domain->giveElement(i)->giveDofManagerNumber(j);
3437  }
3438 
3439  _element = new Subdivision :: RS_Triangle(i, mesh, i, enodes);
3440  } else if ( domain->giveElement(i)->giveGeometryType() == EGT_tetra_1 ) {
3441  enodes.resize(4);
3442  for ( int j = 1; j <= 4; j++ ) {
3443  enodes.at(j) = domain->giveElement(i)->giveDofManagerNumber(j);
3444  }
3445 
3446  _element = new Subdivision :: RS_Tetra(i, mesh, i, enodes);
3447  } else {
3448  OOFEM_ERROR("Unsupported element geometry (element %d)", i);
3449  _element = NULL;
3450  }
3451  _element->setGlobalNumber( domain->giveElement(i)->giveGlobalNumber() );
3452 #ifdef __PARALLEL_MODE
3453  _element->setParallelMode( domain->giveElement(i)->giveParallelMode() );
3454 #endif
3455  this->mesh->addElement(_element);
3456 
3457  }
3458 
3459  // import connectivities for local elements only
3460  for ( int i = 1; i <= nelems; i++ ) {
3461 #ifdef __PARALLEL_MODE
3462  if ( this->mesh->giveElement(i)->giveParallelMode() != Element_local ) {
3463  continue;
3464  }
3465 
3466 #endif
3468  }
3469 
3470 #ifdef __PARALLEL_MODE
3471  // import connectivities for shared nodes only
3472  // build global shared node map (for top level only)
3473  this->mesh->initGlobalSharedNodeMap();
3474  for ( int i = 1; i <= nnodes; i++ ) {
3475  _node = this->mesh->giveNode(i);
3476  if ( _node->giveParallelMode() != DofManager_shared ) {
3477  continue;
3478  }
3479 
3481  this->mesh->insertGlobalSharedNodeMap(_node);
3482  }
3483 
3484  // build array of shared edges (for top level only)
3485  // this can be done after connectivity is available for all shared nodes
3486  for ( int i = 1; i <= nnodes; i++ ) {
3487  _node = this->mesh->giveNode(i);
3488  if ( _node->giveParallelMode() != DofManager_shared ) {
3489  continue;
3490  }
3491 
3492  _node->numberSharedEdges();
3493  }
3494 
3495  // initiate queue for shared edge exchange
3496  for ( int i = 1; i <= mesh->giveNumberOfEdges(); i++ ) {
3497  if ( mesh->giveEdge(i)->givePartitions()->giveSize() ) {
3498  sharedEdgesQueue.push_back(i);
3499  }
3500  }
3501 
3502  // exchange shared edges (on top level)
3503  this->exchangeSharedEdges();
3504 #endif
3505 
3506 #ifdef __OOFEG
3507  #ifdef DRAW_MESH_BEFORE_BISECTION
3508  nelems = mesh->giveNumberOfElements();
3509  for ( int i = 1; i <= nelems; i++ ) {
3510  #ifdef __PARALLEL_MODE
3511  if ( this->mesh->giveElement(i)->giveParallelMode() != Element_local ) {
3512  continue;
3513  }
3514 
3515  #endif
3516  mesh->giveElement(i)->drawGeometry();
3517  }
3518 
3519  ESIEventLoop(YES, "Before bisection; Press Ctrl-p to continue");
3520  #endif
3521 #endif
3522 
3523  // bisect mesh
3524  this->bisectMesh();
3525  // smooth mesh
3526  if ( smoothingFlag ) {
3527  // smooth only if there are new irregulars;
3528  // this is reasonable only if the smoothing is done locally !!!
3529  if ( nnodes != mesh->giveNumberOfNodes() ) {
3530  this->smoothMesh();
3531  }
3532  }
3533 
3534 #ifdef __OOFEG
3535  #ifdef DRAW_MESH_AFTER_BISECTION
3536  nelems = mesh->giveNumberOfElements();
3537  for ( int i = 1; i <= nelems; i++ ) {
3538  if ( !mesh->giveElement(i)->isTerminal() ) {
3539  continue;
3540  }
3541 
3542  #ifdef __PARALLEL_MODE
3543  if ( this->mesh->giveElement(i)->giveParallelMode() != Element_local ) {
3544  continue;
3545  }
3546 
3547  #endif
3548  mesh->giveElement(i)->drawGeometry();
3549  }
3550 
3551  ESIEventLoop(YES, "After bisection; Press Ctrl-p to continue");
3552  #endif
3553 #endif
3554 
3555  Dof *dof;
3556  DofManager *parentNodePtr, *node;
3557  Element *elem;
3558  CrossSection *crossSection;
3559  Material *mat;
3560  NonlocalBarrier *barrier;
3562  InitialCondition *ic;
3563  Function *func;
3564  Set *set;
3565  std :: string name;
3566 
3567  // create new mesh (missing param for new mesh!)
3568  nnodes = mesh->giveNumberOfNodes();
3569  ( * dNew ) = new Domain( 2, domain->giveSerialNumber() + 1, domain->giveEngngModel() );
3570  ( * dNew )->setDomainType( domain->giveDomainType() );
3571 
3572  // copy dof managers
3573  ( * dNew )->resizeDofManagers(nnodes);
3574  const IntArray dofIDArrayPtr = domain->giveDefaultNodeDofIDArry();
3575  for ( int inode = 1; inode <= nnodes; inode++ ) {
3576  parent = mesh->giveNode(inode)->giveParent();
3577  if ( parent ) {
3578  parentNodePtr = domain->giveNode(parent);
3579  // inherit all data from parent (bc, ic, load, etc.)
3580  node = classFactory.createDofManager(parentNodePtr->giveInputRecordName(), inode, * dNew);
3581  node->setNumberOfDofs(0);
3582  node->setLoadArray( * parentNodePtr->giveLoadArray() );
3583  // create individual DOFs
3584  for ( Dof *idofPtr: *parentNodePtr ) {
3585  if ( idofPtr->isPrimaryDof() ) {
3586  dof = new MasterDof( node, idofPtr->giveBcId(), idofPtr->giveIcId(), idofPtr->giveDofID() );
3587  } else {
3588  SimpleSlaveDof *simpleSlaveDofPtr = dynamic_cast< SimpleSlaveDof * >(idofPtr);
3589  if ( simpleSlaveDofPtr ) {
3590  dof = new SimpleSlaveDof( node, simpleSlaveDofPtr->giveMasterDofManagerNum(), idofPtr->giveDofID() );
3591  } else {
3592  OOFEM_ERROR("unsupported DOF type");
3593  dof = NULL;
3594  }
3595  }
3596  node->appendDof(dof);
3597  }
3598 
3599  node->setGlobalNumber( parentNodePtr->giveGlobalNumber() );
3600 #ifdef __PARALLEL_MODE
3601  node->setParallelMode( parentNodePtr->giveParallelMode() );
3602  node->setPartitionList( parentNodePtr->givePartitionList() );
3603 #endif
3604  } else {
3605  // newly created node (irregular)
3606  node = new Node(inode, *dNew);
3607  //create new node with default DOFs
3608  int ndofs = dofIDArrayPtr.giveSize();
3609  node->setNumberOfDofs(0);
3610 
3611  // create individual DOFs
3612  for ( int idof = 1; idof <= ndofs; idof++ ) {
3613 #ifdef NM
3614  dof = NULL;
3615  FloatArray *coords = mesh->giveNode(inode)->giveCoordinates();
3616  if ( !dof ) {
3617  if ( fabs(coords->at(1) - 200.0) < 0.000001 ) {
3618  if ( coords->at(2) > -0.000001 && coords->at(2) < 97.500001 ) {
3619  dof = new MasterDof( node, 1, 0, dofIDArrayPtr.at ( idof ) );
3620  }
3621  }
3622  }
3623 
3624  if ( !dof ) {
3625  if ( fabs( coords->at(2) ) < 0.000001 ) {
3626  dof = new MasterDof( node, 1, 0, dofIDArrayPtr.at ( idof ) );
3627  }
3628  }
3629 
3630  if ( !dof ) {
3631  if ( fabs( coords->at(1) ) < 0.000001 ) {
3632  if ( coords->at(2) > 102.4999999 && coords->at(2) < 199.9999999 ) {
3633  dof = new SimpleSlaveDof( node, 5, ( DofID ) dofIDArrayPtr.at ( idof ) ); // HUHU
3634  }
3635  }
3636  }
3637 
3638  if ( !dof ) {
3639  if ( fabs(coords->at(2) - 200.0) < 0.000001 ) {
3640  if ( coords->at(1) > 0.000001 && coords->at(1) < 200.000001 ) {
3641  dof = new SimpleSlaveDof( node, 6, ( DofID ) dofIDArrayPtr.at ( idof ) ); // HUHU
3642  }
3643  }
3644  }
3645 
3646  if ( !dof ) {
3647  dof = new MasterDof( node, 0, 0, dofIDArrayPtr.at ( idof ) );
3648  }
3649 
3650 #else
3651  #ifdef BRAZIL_2D
3652  dof = NULL;
3653  FloatArray *coords = mesh->giveNode(inode)->giveCoordinates();
3654  if ( !dof ) {
3655  if ( fabs( coords->at(1) ) < 0.000001 ) {
3656  if ( coords->at(2) > 0.000001 ) {
3657  if ( idof == 1 ) {
3658  dof = new MasterDof( node, 1, 0, dofIDArrayPtr.at ( idof ) );
3659  }
3660  }
3661  }
3662  }
3663 
3664  if ( !dof ) {
3665  if ( fabs( coords->at(2) ) < 0.000001 ) {
3666  if ( coords->at(1) > 0.000001 ) {
3667  if ( idof == 2 ) {
3668  dof = new MasterDof( node, 1, 0, dofIDArrayPtr.at ( idof ) );
3669  }
3670  }
3671  }
3672  }
3673 
3674  if ( !dof ) {
3675  dof = new MasterDof( node, 0, 0, dofIDArrayPtr.at ( idof ) );
3676  }
3677 
3678  #else
3679  #ifdef THREEPBT_3D
3680  dof = NULL;
3681  FloatArray *coords = mesh->giveNode(inode)->giveCoordinates();
3682  if ( !dof ) {
3683  if ( fabs( coords->at(1) ) < 0.000001 ) {
3684  if ( coords->at(2) > 0.000001 ) {
3685  if ( coords->at(3) > 499.999999 ) {
3686  if ( idof == 1 || idof == 3 ) {
3687  dof = new MasterDof( node, 1, 0, dofIDArrayPtr.at ( idof ) );
3688  }
3689  }
3690  }
3691  }
3692  }
3693 
3694  if ( !dof ) {
3695  if ( coords->at(1) > 1999.999999 ) {
3696  if ( coords->at(2) > 0.000001 ) {
3697  if ( coords->at(3) > 499.999999 ) {
3698  if ( idof == 3 ) {
3699  dof = new MasterDof( node, 1, 0, dofIDArrayPtr.at ( idof ) );
3700  }
3701  }
3702  }
3703  }
3704  }
3705 
3706  if ( !dof ) {
3707  dof = new MasterDof( node, 0, 0, dofIDArrayPtr.at ( idof ) );
3708  }
3709 
3710  #else
3711  #ifdef HEADEDSTUD
3712  double dist, rad;
3713  FloatArray *coords = mesh->giveNode(inode)->giveCoordinates();
3714 
3715  dof = NULL;
3716  if ( !dof ) {
3717  dist = coords->at(1) * coords->at(1) + coords->at(3) * coords->at(3);
3718  if ( coords->at(2) < 88.000001 ) {
3719  if ( coords->at(2) > 70.000001 ) {
3720  if ( fabs(dist - 7.0 * 7.0) < 0.01 ) { // be very tolerant (geometry is not precise)
3721  if ( idof == 1 || idof == 3 ) {
3722  dof = new MasterDof( node, 1, 0, ( DofIDItem ) dofIDArrayPtr.at ( idof ) );
3723  }
3724  }
3725  } else if ( coords->at(2) > 67.9999999 ) {
3726  rad = 18.0 - 11.0 / 5.5 * ( coords->at(2) - 64.5 );
3727  if ( fabs(dist - rad * rad) < 0.01 ) { // be very tolerant (geometry is not precise)
3728  if ( idof == 1 || idof == 3 ) {
3729  dof = new MasterDof( node, 1, 0, ( DofIDItem ) dofIDArrayPtr.at ( idof ) );
3730  }
3731 
3732  if ( idof == 2 ) {
3733  dof = new MasterDof( node, 2, 0, ( DofIDItem ) dofIDArrayPtr.at ( idof ) );
3734  OOFEM_LOG_INFO("Subdivision: Conrolled Node %d", inode);
3735  }
3736  }
3737  }
3738  }
3739  }
3740 
3741  if ( !dof ) {
3742  if ( coords->at(1) > 299.999999 || coords->at(3) > 299.999999 ) {
3743  dof = new MasterDof( node, 1, 0, ( DofIDItem ) dofIDArrayPtr.at ( idof ) );
3744  }
3745  }
3746 
3747  if ( !dof ) {
3748  if ( coords->at(1) < 0.00000001 ) {
3749  if ( idof == 1 ) {
3750  dof = new MasterDof( node, 1, 0, ( DofIDItem ) dofIDArrayPtr.at ( idof ) );
3751  }
3752  }
3753  }
3754 
3755  if ( !dof ) {
3756  dof = new MasterDof( node, 0, 0, ( DofIDItem ) dofIDArrayPtr.at ( idof ) );
3757  }
3758 
3759  #else
3760  dof = new MasterDof( node, 0, 0, dofIDArrayPtr.at ( idof ) );
3761  #endif
3762  #endif
3763  #endif
3764 #endif
3765  node->appendDof(dof);
3766  }
3767 
3768  node->setGlobalNumber( mesh->giveNode(inode)->giveGlobalNumber() );
3769 #ifdef __PARALLEL_MODE
3770  node->setParallelMode( mesh->giveNode(inode)->giveParallelMode() );
3771  node->setPartitionList( mesh->giveNode(inode)->givePartitions() );
3772 #endif
3773  }
3774 
3775  // set node coordinates
3776  static_cast< Node * >(node)->setCoordinates( * mesh->giveNode(inode)->giveCoordinates() );
3777  node->setBoundaryFlag( mesh->giveNode(inode)->isBoundary() );
3778  ( * dNew )->setDofManager(inode, node);
3779  } // end creating dof managers
3780 
3781  // create elements
3782  // count number of local terminal elements first
3783  int nterminals = 0;
3784  nelems = mesh->giveNumberOfElements();
3785  for ( int ielem = 1; ielem <= nelems; ielem++ ) {
3786 #ifdef __PARALLEL_MODE
3787  if ( mesh->giveElement(ielem)->giveParallelMode() != Element_local ) {
3788  continue;
3789  }
3790 
3791 #endif
3792  if ( mesh->giveElement(ielem)->isTerminal() ) {
3793  nterminals++;
3794  }
3795  }
3796 
3797 #ifdef __PARALLEL_MODE
3798  IntArray parentElemMap(nterminals);
3799 #endif
3800  ( * dNew )->resizeElements(nterminals);
3801  int eNum = 0;
3802  for ( int ielem = 1; ielem <= nelems; ielem++ ) {
3803 #ifdef __PARALLEL_MODE
3804  if ( mesh->giveElement(ielem)->giveParallelMode() != Element_local ) {
3805  continue;
3806  }
3807 
3808 #endif
3809  if ( !mesh->giveElement(ielem)->isTerminal() ) {
3810  continue;
3811  }
3812 
3813  eNum++;
3814  parent = mesh->giveElement(ielem)->giveTopParent();
3815 #ifdef __PARALLEL_MODE
3816  parentElemMap.at(eNum) = parent;
3817 #endif
3818  if ( parent ) {
3819  // Copy most of the existing parent element:
3820  DynamicInputRecord ir( *domain->giveElement ( parent ) );
3822  ir.giveRecordKeywordField(name);
3823  elem = classFactory.createElement(name.c_str(), eNum, * dNew);
3824  elem->initializeFrom(& ir);
3825  elem->setGlobalNumber( mesh->giveElement(ielem)->giveGlobalNumber() );
3826 #ifdef __PARALLEL_MODE
3827  //ir.setRecordKeywordNumber( mesh->giveElement(ielem)->giveGlobalNumber() );
3828  // not subdivided elements inherit globNum, subdivided give -1
3829  // local elements have array partitions empty !
3830 #endif
3831  ( * dNew )->setElement(eNum, elem);
3832  } else {
3833  OOFEM_ERROR("parent element missing");
3834  }
3835  } // end loop over elements
3836 
3837  // create the rest of the model description (BCs, CrossSections, Materials, etc)
3838  // cross sections
3839  int ncrosssect = domain->giveNumberOfCrossSectionModels();
3840  ( * dNew )->resizeCrossSectionModels(ncrosssect);
3841  for ( int i = 1; i <= ncrosssect; i++ ) {
3843  ir.giveRecordKeywordField(name);
3844 
3845  crossSection = classFactory.createCrossSection(name.c_str(), i, * dNew);
3846  crossSection->initializeFrom(& ir);
3847  ( * dNew )->setCrossSection(i, crossSection);
3848  }
3849 
3850  // materials
3851  int nmat = domain->giveNumberOfMaterialModels();
3852  ( * dNew )->resizeMaterials(nmat);
3853  for ( int i = 1; i <= nmat; i++ ) {
3854  DynamicInputRecord ir( *domain->giveMaterial ( i ) );
3855  ir.giveRecordKeywordField(name);
3856 
3857  mat = classFactory.createMaterial(name.c_str(), i, * dNew);
3858  mat->initializeFrom(& ir);
3859  ( * dNew )->setMaterial(i, mat);
3860  }
3861 
3862  // barriers
3863  int nbarriers = domain->giveNumberOfNonlocalBarriers();
3864  ( * dNew )->resizeNonlocalBarriers(nbarriers);
3865  for ( int i = 1; i <= nbarriers; i++ ) {
3867  ir.giveRecordKeywordField(name);
3868 
3869  barrier = classFactory.createNonlocalBarrier(name.c_str(), i, * dNew);
3870  barrier->initializeFrom(& ir);
3871  ( * dNew )->setNonlocalBarrier(i, barrier);
3872  }
3873 
3874  // boundary conditions
3876  ( * dNew )->resizeBoundaryConditions(nbc);
3877  for ( int i = 1; i <= nbc; i++ ) {
3878  DynamicInputRecord ir( *domain->giveBc ( i ) );
3879  ir.giveRecordKeywordField(name);
3880 
3881 
3882  bc = classFactory.createBoundaryCondition(name.c_str(), i, * dNew);
3883  bc->initializeFrom(& ir);
3884  ( * dNew )->setBoundaryCondition(i, bc);
3885  }
3886 
3887  // initial conditions
3889  ( * dNew )->resizeInitialConditions(nic);
3890  for ( int i = 1; i <= nic; i++ ) {
3891  DynamicInputRecord ir( *domain->giveIc ( i ) );
3892 
3893  ic = new InitialCondition(i, *dNew);
3894  ic->initializeFrom(& ir);
3895  ( * dNew )->setInitialCondition(i, ic);
3896  }
3897 
3898  // load time functions
3899  int nfunc = domain->giveNumberOfFunctions();
3900  ( * dNew )->resizeFunctions(nfunc);
3901  for ( int i = 1; i <= nfunc; i++ ) {
3902  DynamicInputRecord ir( *domain->giveFunction ( i ) );
3903  ir.giveRecordKeywordField(name);
3904 
3905  func = classFactory.createFunction(name.c_str(), i, * dNew);
3906  func->initializeFrom(& ir);
3907  ( * dNew )->setFunction(i, func);
3908  }
3909 
3910  // sets
3911  int nset = domain->giveNumberOfSets();
3912  ( * dNew )->resizeSets(nset);
3913  for ( int i = 1; i <= nset; i++ ) {
3914  DynamicInputRecord ir( *domain->giveSet ( i ) );
3915  ir.giveRecordKeywordField(name);
3916 
3917  set = new Set(i, * dNew);
3918  set->initializeFrom(& ir);
3919  ( * dNew )->setSet(i, set);
3920  }
3921 
3922 
3923  // post initialize components
3924  ( * dNew )->postInitialize();
3925 
3926  // copy output manager settings
3927  ( * dNew )->giveOutputManager()->beCopyOf( domain->giveOutputManager() );
3928 
3929  timer.stopTimer();
3930 #ifdef __PARALLEL_MODE
3931  OOFEM_LOG_INFO( "[%d] Subdivision: created new mesh (%d nodes and %d elements) in %.2fs\n",
3932  ( * dNew )->giveEngngModel()->giveRank(), nnodes, eNum, timer.getUtime() );
3933 #else
3934  OOFEM_LOG_INFO( "Subdivision: created new mesh (%d nodes and %d elements) in %.2fs\n",
3935  nnodes, eNum, timer.getUtime() );
3936 #endif
3937 
3938 
3939 
3940 
3941 #ifdef __PARALLEL_MODE
3942  #ifdef __VERBOSE_PARALLEL
3943  nnodes = ( * dNew )->giveNumberOfDofManagers();
3944  for ( int inode = 1; inode <= nnodes; inode++ ) {
3945  if ( ( * dNew )->giveDofManager(inode)->giveParallelMode() == DofManager_shared ) {
3946  //OOFEM_LOG_INFO ("[%d] Shared Node %d[%d]\n", this->giveRank(), inode, (*dNew)->giveDofManager(inode)->giveGlobalNumber());
3947  }
3948  }
3949 
3950  #endif
3951 
3952  // we need to assign global numbers to newly generated elements
3953  this->assignGlobalNumbersToElements(* dNew);
3954 
3955  bool nonloc = false;
3956  nmat = ( * dNew )->giveNumberOfMaterialModels();
3957  for ( int im = 1; im <= nmat; im++ ) {
3958  if ( ( * dNew )->giveMaterial(im)->giveInterface(NonlocalMaterialExtensionInterfaceType) ) {
3959  nonloc = true;
3960  }
3961  }
3962 
3963  if ( nonloc ) {
3964  exchangeRemoteElements(* dNew, parentElemMap);
3965  }
3966 
3967  ( * dNew )->commitTransactions( ( * dNew )->giveTransactionManager() );
3968 
3969  // print some statistics
3970  if (0) {
3971  for (int in=1; in<=(*dNew)->giveNumberOfDofManagers(); in++) {
3972  DynamicInputRecord ir;
3973  (*dNew)->giveDofManager(in)->giveInputRecord(ir);
3974  OOFEM_LOG_INFO("[%d]:[%d]:%s\n", this->giveRank(), (*dNew)->giveDofManager(in)->giveGlobalNumber(), ir.giveRecordAsString().c_str());
3975  }
3976  IntArray nodes;
3977  for (int in=1; in<=(*dNew)->giveNumberOfElements(); in++) {
3978  DynamicInputRecord ir;
3979  (*dNew)->giveElement(in)->giveInputRecord(ir);
3980  nodes = (*dNew)->giveElement(in)->giveDofManArray();
3981  // translate local node numbers to globnums
3982  for (int ii=1; ii<=nodes.giveSize(); ii++)
3983  nodes.at(ii) = (*dNew)->giveNode(nodes.at(ii))->giveGlobalNumber();
3984  ir.setField(nodes, "gnodes");
3985  OOFEM_LOG_INFO("[%d]:[%d]:%s\n", this->giveRank(), (*dNew)->giveElement(in)->giveGlobalNumber(), ir.giveRecordAsString().c_str());
3986  }
3987  }
3988  nelems = ( * dNew )->giveNumberOfElements();
3989  int localVals [ 2 ], globalVals [ 2 ];
3990  localVals [ 0 ] = 0;
3991  for ( int ielem = 1; ielem <= nelems; ielem++ ) {
3992  if ( ( * dNew )->giveElement(ielem)->giveParallelMode() == Element_local ) {
3993  localVals [ 0 ]++;
3994  }
3995  }
3996 
3997  nnodes = ( * dNew )->giveNumberOfDofManagers();
3998  localVals [ 1 ] = 0;
3999  for ( int inode = 1; inode <= nnodes; inode++ ) {
4000  if ( ( * dNew )->giveDofManager(inode)->isLocal() ) {
4001  localVals [ 1 ]++;
4002  }
4003  }
4004 
4005  MPI_Reduce(localVals, globalVals, 2, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
4006  if ( this->giveRank() == 0 ) {
4007  OOFEM_LOG_INFO("Subdivision: new mesh info: %d nodes, %d elements in total\n", globalVals [ 1 ], globalVals [ 0 ]);
4008  }
4009 #else
4010  // we need to assign global numbers to newly generated elements
4011  this->assignGlobalNumbersToElements(* dNew);
4012 #endif
4013 
4014  return MI_OK;
4015 }
4016 
4017 
4018 void
4020 {
4021  int ie, nelems = mesh->giveNumberOfElements(), nelems_old = 0, terminal_local_elems = nelems;
4022  int nnodes = mesh->giveNumberOfNodes(), nnodes_old;
4023  double iedensity, rdensity;
4024  int repeat = 1, loop = 0, max_loop = 0; // max_loop != 0 use only for debugging
4025  RS_Element *elem;
4026  RS_Node *node;
4027  //std::queue<int> subdivqueue;
4028 #ifdef __PARALLEL_MODE
4029  int remote_elems = 0;
4030  int myrank = this->giveRank();
4031  int problem_size = this->giveNumberOfProcesses();
4032  int value;
4033  int *partitionsIrregulars = new int[ problem_size ];
4034 #endif
4035 
4036  // get the max globnum on the mesh
4037  // determine max global number of local nodes
4038  int maxlocalglobal = 0, maxglobalnumber;
4039  for ( int in = 1; in <= nnodes; in++ ) {
4040  maxlocalglobal = max( maxlocalglobal, mesh->giveNode(in)->giveGlobalNumber() );
4041  }
4042 #ifdef __PARALLEL_MODE
4043  // determine max global number on all partitions
4044  MPI_Allreduce(& maxlocalglobal, & maxglobalnumber, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
4045 #else
4046  maxglobalnumber=maxlocalglobal;
4047 #endif
4048 
4049 
4050  // repeat bisection until no new element is created
4051  while ( repeat && ( loop < max_loop || max_loop == 0 ) ) {
4052 
4053 
4054  nnodes_old = nnodes;
4055 #ifdef __PARALLEL_MODE
4056  OOFEM_LOG_INFO("[%d] Subdivision::bisectMesh: entering bisection loop %d\n", myrank, ++loop);
4057 #else
4058  OOFEM_LOG_INFO("Subdivision::bisectMesh: entering bisection loop %d\n", ++loop);
4059 #endif
4060  repeat = 0;
4061  // process only newly created elements in pass 2 and more
4062  for ( ie = nelems_old + 1; ie <= nelems; ie++ ) {
4063  elem = mesh->giveElement(ie);
4064 #ifdef __PARALLEL_MODE
4065  // skip bisection of remote elements (in first pass only (nelems_old = 0));
4066  // in pass 2 and more there should be no remote elements because
4067  // only newly created elements are processed
4068  if ( nelems_old == 0 ) {
4069  if ( elem->giveParallelMode() == Element_remote ) {
4070  remote_elems++;
4071  terminal_local_elems--;
4072  continue;
4073  }
4074  }
4075 
4076 #endif
4077  if ( !elem->isTerminal() ) {
4078  continue;
4079  }
4080 
4081 #ifdef __PARALLEL_MODE
4082  // bisect local elements only
4083  if ( elem->giveParallelMode() != Element_local ) {
4084  continue;
4085  }
4086 
4087 #endif
4088 
4089 #ifdef DEBUG_CHECK
4090  if ( elem->giveQueueFlag() ) {
4091  OOFEM_ERROR("unexpected queue flag of %d ", ie);
4092  }
4093 
4094 #endif
4095 
4096  iedensity = elem->giveDensity();
4097  rdensity = elem->giveRequiredDensity();
4098 
4099  // first select all candidates for local bisection based on required mesh density
4100 
4101  if ( rdensity < iedensity ) {
4102  subdivqueue.push(ie);
4103  elem->setQueueFlag(true);
4104 
4105 
4106 
4107 #ifndef __PARALLEL_MODE
4108  repeat = 1; // force repetition in seqeuntial run
4109 #endif
4110  /*
4111  * #ifdef __PARALLEL_MODE
4112  * OOFEM_LOG_INFO("[%d] Subdivision: scheduling element %d[%d] for bisection, dens=%lf rdens=%lf\n", myrank, ie, elem->giveGlobalNumber(), iedensity, rdensity);
4113  ******#else
4114  * OOFEM_LOG_INFO("Subdivision: scheduling element %d for bisection, dens=%lf rdens=%lf\n", ie, iedensity, rdensity);
4115  ******#endif
4116  */
4117  }
4118  }
4119 
4120 #ifdef DEBUG_INFO
4121  #ifdef __PARALLEL_MODE
4122  OOFEM_LOG_INFO("[%d] (with %d nodes and %d elems)\n", myrank, nnodes_old, terminal_local_elems + remote_elems);
4123  #else
4124  OOFEM_LOG_INFO("(with %d nodes and %d elems)\n", nnodes_old, terminal_local_elems);
4125  #endif
4126 #endif
4127 
4128 #ifdef __PARALLEL_MODE
4129  for ( value = 0; value == 0; value = exchangeSharedIrregulars() ) {
4130 #endif
4131  // loop over subdivision queue to bisect all local elements there
4132  while ( !subdivqueue.empty() ) {
4133  elem = mesh->giveElement( subdivqueue.front() );
4134 #ifdef DEBUG_CHECK
4135  #ifdef __PARALLEL_MODE
4136  if ( elem->giveParallelMode() != Element_local ) {
4137  OOFEM_ERROR( "nonlocal element %d not expected for bisection", elem->giveNumber() );
4138  }
4139 
4140  #endif
4141 #endif
4142  elem->evaluateLongestEdge();
4144  subdivqueue.pop();
4145  }
4146 
4147 #ifdef __PARALLEL_MODE
4148  // in parallel communicate with neighbours the irregular nodes on shared bondary
4149  }
4150 
4151 #endif
4152 
4153  int in;
4154  nnodes = mesh->giveNumberOfNodes();
4155 
4156 #ifndef __PARALLEL_MODE
4157  int myrank = 0;
4158 #endif
4159  // assign global numbers to newly introduced irregulars while
4160  // keeping global numbering of existing (master) nodes
4161  // idea: first determine the max globnum already assigned
4162  // and start global numbering of new nodes from this value up
4163 
4164  // determine number of local irregulars
4165  // this is needed to determine offsets on each partition to start global numbering
4166  // the numbering of irregulars starts on rank 0 partition by numbering subseqeuntly all its irregulars
4167  // then we proceed with rank 1, etc.
4168  // the shared irregulars receive their number from partition with the lowest rank.
4169 
4170  // count local irregulars that receive their global number from this partition
4171  int localIrregulars = 0;
4172  for ( in = nnodes_old+1; in <= nnodes; in++ ) {
4173  if ( this->isNodeLocalIrregular(mesh->giveNode(in), myrank) && ( mesh->giveNode(in)->giveGlobalNumber() == 0 ) ) {
4174  localIrregulars++;
4175  }
4176  }
4177 
4178 #ifdef __PARALLEL_MODE
4179  #ifdef __VERBOSE_PARALLEL
4180  OOFEM_LOG_INFO("[%d] Subdivision::bisectMesh: number of new local irregulars is %d\n", myrank, localIrregulars);
4181  #endif
4182  int localOffset = 0;
4183  int irank, gnum, globalIrregulars = 0;
4184  // gather number of local irregulars from all partitions
4185  MPI_Allgather(& localIrregulars, 1, MPI_INT, partitionsIrregulars, 1, MPI_INT, MPI_COMM_WORLD);
4186  // compute local offset
4187  for ( irank = 0; irank < myrank; irank++ ) {
4188  localOffset += partitionsIrregulars [ irank ];
4189  }
4190  // start to assign global numbers to local irregulars
4191  int availGlobNum = maxglobalnumber + localOffset;
4192  for ( in = nnodes_old+1; in <= nnodes; in++ ) {
4193  node = mesh->giveNode(in);
4194  if ( this->isNodeLocalIrregular(node, myrank) && ( node->giveGlobalNumber() == 0 ) ) {
4195  // set negative globnum to mark newly assigned nodes with globnum to participate in shared globnum data exchange
4196  node->setGlobalNumber( -( ++availGlobNum ) );
4197  #ifdef __VERBOSE_PARALLEL
4198  //OOFEM_LOG_INFO ("[%d] %d --> [%d]\n", myrank, in, availGlobNum);
4199  #endif
4200  }
4201  }
4202 #else
4203  int localOffset=0;
4204  // start to assign global numbers to local irregulars
4205  int availGlobNum = maxglobalnumber+localOffset;
4206  for ( in = nnodes_old+1; in <= nnodes; in++ ) {
4207  node = mesh->giveNode(in);
4208  node->setGlobalNumber( ( ++availGlobNum ) );
4209  }
4210 
4211 #endif
4212 
4213 #ifdef __PARALLEL_MODE
4214  // finally, communicate global numbers assigned to shared irregulars
4216  for ( in = nnodes_old+1; in <= nnodes; in++ ) {
4217  node = mesh->giveNode(in);
4218  gnum = node->giveGlobalNumber();
4219  if ( gnum < 0 ) {
4220  // turn all globnums to positive
4221  node->setGlobalNumber(-gnum);
4222  // update global shared node map (for refined level)
4223  // this must be done after globnum is made positive
4224  this->mesh->insertGlobalSharedNodeMap(node);
4225  } else if ( gnum == 0 ) {
4226  OOFEM_ERROR("rank [%d] zero globnum identified on node %d", myrank, in);
4227  }
4228  }
4229 
4230  // update max globnum
4231  for ( irank = 0; irank < problem_size; irank++ ) {
4232  globalIrregulars += partitionsIrregulars [ irank ];
4233  }
4234 
4235  if ( globalIrregulars ) {
4236  maxglobalnumber += globalIrregulars;
4237  // exchange shared edges
4238  // this must be done after globnums are assigned to new shared irregulars
4239  /* BP
4240  this->exchangeSharedEdges();
4241  repeat = 1; // force repetition in parallel run
4242  */
4243  }
4244 
4245 #else
4246  maxglobalnumber+=localIrregulars;
4247 #endif
4248 
4249  // symbolic bisection is finished;
4250  // now we need to create new mesh;
4251  // also there is a need to update element connectivities
4252  for ( ie = 1; ie <= nelems; ie++ ) {
4253  elem = mesh->giveElement(ie);
4254  if ( !elem->isTerminal() ) {
4255  continue;
4256  }
4257 
4258 #ifdef __PARALLEL_MODE
4259  if ( elem->giveParallelMode() != Element_local ) {
4260  continue;
4261  }
4262 
4263 #endif
4264  elem->generate(sharedEdgesQueue);
4265  }
4266 
4267 #ifdef __PARALLEL_MODE
4268  if ( globalIrregulars ) {
4269  // exchange shared edges
4270  // this must be done after globnums are assigned to new shared irregulars
4271  this->exchangeSharedEdges();
4272  repeat = 1; // force repetition in parallel run
4273  }
4274 #endif
4275  // unmark local unshared irregulars marked for connectivity setup
4276  // important this may be not done before generate !!!
4277  for ( in = nnodes_old+1; in <= nnodes; in++ ) {
4278  if ( mesh->giveNode(in)->giveNumber() < 0 ) {
4279  mesh->giveNode(in)->setNumber( -mesh->giveNode(in)->giveNumber() );
4280  }
4281  }
4282 
4283  nelems_old = nelems;
4284  nelems = mesh->giveNumberOfElements();
4285  terminal_local_elems = 0;
4286  for ( ie = 1; ie <= nelems; ie++ ) {
4287  elem = mesh->giveElement(ie);
4288  if ( !elem->isTerminal() ) {
4289  continue;
4290  }
4291 
4292 #ifdef __PARALLEL_MODE
4293  if ( elem->giveParallelMode() != Element_local ) {
4294  continue;
4295  }
4296 
4297 #endif
4298  elem->update_neighbours();
4299  terminal_local_elems++;
4300  }
4301 
4302 #if 0
4303  #ifdef __PARALLEL_MODE
4304  int global_repeat = 0;
4305 
4306  // determine whether any partition needs additional bisection pass
4307  MPI_Allreduce(& repeat, & global_repeat, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
4308  repeat = global_repeat;
4309  #endif
4310 #endif
4311 
4312 #ifdef __OOFEG
4313  #ifdef DRAW_MESH_AFTER_EACH_BISECTION_LEVEL
4314  for ( ie = 1; ie <= nelems; ie++ ) {
4315  if ( !mesh->giveElement(ie)->isTerminal() ) {
4316  continue;
4317  }
4318 
4319  mesh->giveElement(ie)->drawGeometry();
4320  }
4321 
4322  ESIEventLoop( YES, const_cast< char * >("Subdivision Bisection; Press Ctrl-p to continue") );
4323  #endif
4324 #endif
4325  } // end bisection loop
4326 #ifdef __PARALLEL_MODE
4327  if (partitionsIrregulars) {
4328  delete[] partitionsIrregulars;
4329  }
4330 #endif
4331 }
4332 
4333 
4334 void
4336 {
4337  int nnodes, nelems, i, j, in, ie;
4338  int pos, number, reg, nd, cycles = 6;
4339  IntArray snodes;
4340  FloatArray *coords;
4341  RS_Element *elem;
4342  //bool fixed;
4343  IntArray node_num_elems, node_con_elems;
4344  IntArray node_num_nodes, node_con_nodes;
4345 
4346 #ifdef DEBUG_SMOOTHING
4347  char buffer [ 1024 ];
4348  int len;
4349 #endif
4350 #ifdef __PARALLEL_MODE
4351  OOFEM_LOG_INFO( "[%d] Subdivision::smoothMesh\n", this->giveRank() );
4352 #else
4353  OOFEM_LOG_INFO("Subdivision::smoothMesh\n");
4354 #endif
4355  nnodes = mesh->giveNumberOfNodes();
4356  nelems = mesh->giveNumberOfElements();
4357 
4358  // count number of elements incident to nodes
4359  node_num_elems.resize(nnodes + 1);
4360  node_num_elems.zero();
4361 
4362  for ( ie = 1; ie <= nelems; ie++ ) {
4363  elem = mesh->giveElement(ie);
4364  if ( !elem->isTerminal() ) {
4365  continue;
4366  }
4367 
4368 #ifdef __PARALLEL_MODE
4369  if ( elem->giveParallelMode() != Element_local ) {
4370  continue;
4371  }
4372 
4373 #endif
4374 
4375  for ( i = 1; i <= elem->giveNodes()->giveSize(); i++ ) {
4376  node_num_elems.at( elem->giveNode(i) )++;
4377  }
4378  }
4379 
4380  // recalculate the number of elements to current addresses
4381  pos = 1;
4382  for ( in = 1; in <= nnodes; in++ ) {
4383  number = node_num_elems.at(in);
4384  node_num_elems.at(in) = pos;
4385  pos += number;
4386  }
4387 
4388  node_num_elems.at(nnodes + 1) = pos;
4389  node_con_elems.resize(pos);
4390 
4391  // store element numbers incident to nodes
4392  for ( ie = 1; ie <= nelems; ie++ ) {
4393  elem = mesh->giveElement(ie);
4394  if ( !elem->isTerminal() ) {
4395  continue;
4396  }
4397 
4398 #ifdef __PARALLEL_MODE
4399  if ( elem->giveParallelMode() != Element_local ) {
4400  continue;
4401  }
4402 
4403 #endif
4404 
4405  for ( i = 1; i <= elem->giveNodes()->giveSize(); i++ ) {
4406  node_con_elems.at(node_num_elems.at( elem->giveNode(i) )++) = ie;
4407  }
4408  }
4409 
4410  // recalculate the addresses to address of the first element
4411  pos = 1;
4412  for ( in = 1; in <= nnodes; in++ ) {
4413  number = node_num_elems.at(in) - pos;
4414  node_num_elems.at(in) = pos;
4415  pos += number;
4416  }
4417 
4418 #ifdef DEBUG_SMOOTHING
4419  #ifdef __PARALLEL_MODE
4420  OOFEM_LOG_INFO( "[%d] Subdivision::smoothMesh: connectivity node_id: elem_ids\n", this->giveRank() );
4421  #else
4422  OOFEM_LOG_INFO("Subdivision::smoothMesh: connectivity node_id: elem_ids\n");
4423  #endif
4424  for ( in = 1; in <= nnodes; in++ ) {
4425  len = 0;
4426  for ( i = node_num_elems.at(in); i < node_num_elems.at(in + 1); i++ ) {
4427  sprintf( buffer + len, " %d", node_con_elems.at(i) );
4428  len = strlen(buffer);
4429  if ( len >= 1024 ) {
4430  OOFEM_ERROR("too small buffer");
4431  }
4432  }
4433 
4434  OOFEM_LOG_INFO("%d:%s\n", in, buffer);
4435  }
4436 
4437 #endif
4438 
4439  // count number of nodes incident to nodes
4440  node_num_nodes.resize(nnodes + 1);
4441  node_num_nodes.zero();
4442 
4443  for ( in = 1; in <= nnodes; in++ ) {
4444  for ( i = node_num_elems.at(in); i < node_num_elems.at(in + 1); i++ ) {
4445  elem = mesh->giveElement( node_con_elems.at(i) );
4446  for ( j = 1; j <= elem->giveNodes()->giveSize(); j++ ) {
4447  nd = elem->giveNode(j);
4448  if ( nd != in && mesh->giveNode(nd)->giveNumber() > 0 ) {
4449  mesh->giveNode(nd)->setNumber(-nd); // mark connected node
4450  node_num_nodes.at(in)++;
4451  }
4452  }
4453  }
4454 
4455  // unmark connected nodes
4456  for ( i = node_num_elems.at(in); i < node_num_elems.at(in + 1); i++ ) {
4457  elem = mesh->giveElement( node_con_elems.at(i) );
4458  for ( j = 1; j <= elem->giveNodes()->giveSize(); j++ ) {
4459  nd = elem->giveNode(j);
4460  mesh->giveNode(nd)->setNumber(nd);
4461  }
4462  }
4463  }
4464 
4465  // recalculate the number of nodes to current addresses
4466  pos = 1;
4467  for ( in = 1; in <= nnodes; in++ ) {
4468  number = node_num_nodes.at(in);
4469  node_num_nodes.at(in) = pos;
4470  pos += number;
4471  }
4472 
4473  node_num_nodes.at(nnodes + 1) = pos;
4474  node_con_nodes.resize(pos);
4475 
4476  // store nodes incident to nodes
4477  for ( in = 1; in <= nnodes; in++ ) {
4478  for ( i = node_num_elems.at(in); i < node_num_elems.at(in + 1); i++ ) {
4479  elem = mesh->giveElement( node_con_elems.at(i) );
4480  for ( j = 1; j <= elem->giveNodes()->giveSize(); j++ ) {
4481  nd = elem->giveNode(j);
4482  if ( nd != in && mesh->giveNode(nd)->giveNumber() > 0 ) {
4483  mesh->giveNode(nd)->setNumber(-nd); // mark connected node
4484  node_con_nodes.at(node_num_nodes.at(in)++) = nd;
4485  }
4486  }
4487  }
4488 
4489  // unmark connected nodes
4490  for ( i = node_num_elems.at(in); i < node_num_elems.at(in + 1); i++ ) {
4491  elem = mesh->giveElement( node_con_elems.at(i) );
4492  for ( j = 1; j <= elem->giveNodes()->giveSize(); j++ ) {
4493  nd = elem->giveNode(j);
4494  mesh->giveNode(nd)->setNumber(nd);
4495  }
4496  }
4497  }
4498 
4499  // recalculate the addresses to address of the first node
4500  pos = 1;
4501  for ( in = 1; in <= nnodes; in++ ) {
4502  number = node_num_nodes.at(in) - pos;
4503  node_num_nodes.at(in) = pos;
4504  pos += number;
4505  }
4506 
4507 #ifdef DEBUG_SMOOTHING
4508  #ifdef __PARALLEL_MODE
4509  OOFEM_LOG_INFO( "[%d] Subdivision::smoothMesh: connectivity node_id: node_ids\n", this->giveRank() );
4510  #else
4511  OOFEM_LOG_INFO("Subdivision::smoothMesh: connectivity node_id: node_ids\n");
4512  #endif
4513  for ( in = 1; in <= nnodes; in++ ) {
4514  len = 0;
4515  for ( i = node_num_nodes.at(in); i < node_num_nodes.at(in + 1); i++ ) {
4516  sprintf( buffer + len, " %d", node_con_nodes.at(i) );
4517  len = strlen(buffer);
4518  if ( len >= 1024 ) {
4519  OOFEM_ERROR("too small buffer");
4520  }
4521  }
4522 
4523  OOFEM_LOG_INFO("%d:%s\n", in, buffer);
4524  }
4525 
4526 #endif
4527 
4528  // identify fixed nodes which should not be subjected to smoothing;
4529  // these are nodes on boundary (geometrical or material);
4530  // note that some of these node may be already marked as boundary if this flag was available on input
4531  // or when processed during bisection or when processed (and assigned) during previous smoothing;
4532  // the last case is applicable only if node's method storeYourself packs the boundary flag;
4533  bool fixed;
4534  for ( in = 1; in <= nnodes; in++ ) {
4535  if ( mesh->giveNode(in)->isBoundary() ) {
4536  continue; // skip boundary node
4537  }
4538 
4539  if ( node_num_elems.at(in) == node_num_elems.at(in + 1) ) {
4540  continue; // skip isolated node
4541  }
4542 
4543  fixed = false;
4544  elem = mesh->giveElement( node_con_elems.at( node_num_elems.at(in) ) );
4545  // check for geometrical boundary
4546  for ( j = 1; j <= elem->giveNeighbors()->giveSize(); j++ ) {
4547  if ( elem->giveNeighbor(j) == 0 ) {
4548  elem->giveSideNodes(j, snodes);
4549  if ( snodes.findFirstIndexOf(in) != 0 ) {
4550  mesh->giveNode(in)->setNumber(-in); // mark fixed node on geometrical boundary
4551  fixed = true;
4552  break;
4553  }
4554  }
4555  }
4556 
4557  if ( !fixed ) {
4558  reg = domain->giveElement( elem->giveTopParent() )->giveRegionNumber();
4559  for ( i = node_num_elems.at(in) + 1; i < node_num_elems.at(in + 1); i++ ) {
4560  elem = mesh->giveElement( node_con_elems.at(i) );
4561  // check for material boundary
4562  if ( domain->giveElement( elem->giveTopParent() )->giveRegionNumber() != reg ) {
4563  mesh->giveNode(in)->setNumber(-in); //mark fixed node on material boundary
4564  fixed = true;
4565  break;
4566  }
4567 
4568  // check for geometrical boundary
4569  for ( j = 1; j <= elem->giveNeighbors()->giveSize(); j++ ) {
4570  if ( elem->giveNeighbor(j) == 0 ) {
4571  elem->giveSideNodes(j, snodes);
4572  if ( snodes.findFirstIndexOf(in) != 0 ) {
4573  mesh->giveNode(in)->setNumber(-in); //mark fixed node on geometrical boundary
4574  fixed = true;
4575  break;
4576  }
4577  }
4578  }
4579 
4580  if ( fixed ) {
4581  break;
4582  }
4583  }
4584  }
4585  }
4586 
4587 #ifdef QUICK_HACK
4588  int jn;
4589  IntArray orderedNodes;
4591  orderedNodes.resize(nnodes);
4592  for ( jn = 1; jn <= nnodes; jn++ ) {
4593  orderedNodes.at(jn) = jn;
4594  }
4595 
4596  sort(orderedNodes, cmp);
4597 #endif
4598 
4599  while ( cycles-- ) {
4600 #ifdef QUICK_HACK
4601  for ( jn = 1; jn <= nnodes; jn++ ) {
4602  in = orderedNodes.at(jn);
4603 #else
4604  for ( in = 1; in <= nnodes; in++ ) {
4605 #endif
4606 #ifdef __PARALLEL_MODE
4607  if ( ( mesh->giveNode(in)->giveParallelMode() == DofManager_shared ) ||
4608  ( mesh->giveNode(in)->giveParallelMode() == DofManager_null ) ) {
4609  continue; // skip shared and remote node
4610  }
4611 
4612 #endif
4613  if ( mesh->giveNode(in)->giveNumber() < 0 ) {
4614  continue; // skip fixed node
4615  }
4616 
4617  if ( mesh->giveNode(in)->isBoundary() ) {
4618  continue; // skip boundary node
4619  }
4620 
4621  coords = mesh->giveNode(in)->giveCoordinates();
4622 
4623 #ifdef DEBUG_CHECK
4624  if ( coords ) {
4625  int count = 0;
4626  coords->zero();
4627  for ( i = node_num_nodes.at(in); i < node_num_nodes.at(in + 1); i++ ) {
4628  if ( mesh->giveNode( node_con_nodes.at(i) ) ) {
4629  if ( mesh->giveNode( node_con_nodes.at(i) )->giveCoordinates() ) {
4630  coords->add( *(mesh->giveNode( node_con_nodes.at(i))->giveCoordinates()));
4631  count++;
4632  } else {
4633  OOFEM_ERROR("node %d without coordinates", in);
4634  }
4635  } else {
4636  OOFEM_ERROR("undefined node %d", in);
4637  }
4638  }
4639 
4640  if ( !count ) {
4641  OOFEM_ERROR("node %d without connectivity", in);
4642  }
4643 
4644  coords->times( 1.0 / ( node_num_nodes.at(in + 1) - node_num_nodes.at(in) ) );
4645  } else {
4646  OOFEM_ERROR("node %d without coordinates", in);
4647  }
4648 
4649 #else
4650  coords->zero();
4651  for ( i = node_num_nodes.at(in); i < node_num_nodes.at(in + 1); i++ ) {
4652  coords->add( * mesh->giveNode( node_con_nodes.at(i) )->giveCoordinates() );
4653  }
4654 
4655  coords->times( 1.0 / ( node_num_nodes.at(in + 1) - node_num_nodes.at(in) ) );
4656 #endif
4657  }
4658  }
4659 
4660  // unmark