MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
schurdataman.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/schurdataman.cc,v 1.84 2017/01/12 14:46:10 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 /* Schur Data Manager */
33 
34 /*
35  * Copyright 1999-2017 Giuseppe Quaranta <quaranta@aero.polimi.it>
36  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
37  */
38 
39 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
40 
41 /* libraries for mesh partitioning computation */
42 #include "metiswrap.h"
43 #include "chacowrap.h"
44 
45 #include "schurdataman.h"
46 #include "mbcomm.h"
47 #include "mysleep.h"
48 #include "except.h"
49 #include "solver.h"
50 #undef min
51 #undef max
52 #include <vector>
53 #include <algorithm>
54 
55 #include "rotor.h"
56 
57 /* struttura contenente le distribuzioni di dofs */
58 struct Adjacency{
59  int *pXadj; /* x ogni nodo: pos. su pAdjncy lista connessioni */
60  int *pAdjncy; /* lista liste connessioni */
61 };
62 
63 /* !!!!!!! Costanti da verificare !!!!!*/
67 const int ADJ_UNDEFINED = -1;
68 
69 #ifndef USE_MPI
71  unsigned OF,
72  Solver* pS,
73  doublereal dInitialTime,
74  const char* sOutputFileName,
75  const char* sInputFileName,
76  bool bAbortAfterInput)
77 : DataManager(HP, OF, pS, dInitialTime, sOutputFileName, sInputFileName, bAbortAfterInput)
78 {
79  silent_cerr("fatal error: you are building SchurDataManager, "
80  "but mbdyn was compiled without MPI. "
81  "Something weird is happening. "
82  "Anyway, please compile with -DUSE_MPI "
83  "to enable parallel solution" << std::endl);
85 }
86 
88 {
89  NO_OP;
90 }
91 
92 void
94 {
95  NO_OP;
96 }
97 
98 void
100 {
101  NO_OP;
102 }
103 
104 void
106 {
107  NO_OP;
108 }
109 
110 void
112  VectorHandler&, VectorHandler&) const
113 {
114  NO_OP;
115 }
116 
117 void
119 {
120  NO_OP;
121 }
122 
123 void
125 {
126  NO_OP;
127 }
128 
129 void
131 {
132  NO_OP;
133 }
134 
135 integer
137 {
138  return 0;
139 }
140 
141 integer *
143 {
144  return 0;
145 }
146 
147 Dof *
149 {
150  return 0;
151 }
152 
153 void
154 SchurDataManager::Output(bool force) const
155 {
156  NO_OP;
157 }
158 #else /* USE_MPI */
159 
160 /* NOTE: define to use Wait/Waitall instead of Test/Testall
161  * Apparently, this results in far better performances,
162  * so we might want to extend it to all other communications */
163 #define USE_MPI_WAIT
164 
165 /* Costruttore - begin */
167  unsigned OF,
168  Solver *pS,
169  doublereal dInitialTime,
170  const char* sOutputFileName,
171  const char* sInputFileName,
172  bool bAbortAfterInput)
173 :DataManager(HP, OF, pS, dInitialTime, sOutputFileName, sInputFileName, bAbortAfterInput),
174 iTotVertices(0),
175 ppMyElems(NULL),
176 iNumLocElems(0),
177 ppMyIntElems(NULL),
178 iNumIntElems(0),
179 ppMyNodes(NULL),
180 iNumLocNodes(0),
181 pLocalDofs(NULL),
182 iNumLocDofs(0),
183 pLocalIntDofs(NULL),
184 iNumIntDofs(0),
185 ppIntNodes(NULL),
186 iNumIntNodes(0),
187 iNumMyInt(0),
188 pMyIntDofs(NULL),
189 pLabelsList(NULL),
190 wgtflag(WEIGHT_VERTICES),
191 pParAmgProcs(NULL),
192 Partitioner(PARTITIONER_DEFAULT),
193 pIndVelComm(NULL),
194 ppExpCntNodes(NULL),
195 ppExpCntElems(NULL),
196 iTotalExpConnections(0)
197 {
198  DEBUGCOUT("Entering SchurDataManager" << std::endl);
199 
200  /* Inizializza il communicator */
201  DataComm = MBDynComm.Dup();
202  DataCommSize = DataComm.Get_size();
203  MyRank = DataComm.Get_rank();
204 
205  DEBUGCOUT("Communicator Size: " << DataCommSize << std::endl);
206 
207  iTotVertices = iTotNodes + Elems.size();
208  DEBUGCOUT("iTotVertices: " << iTotVertices << std::endl);
209 
210  /* parole chiave del blocco parallelizzazione */
211  const char* sKeyWords[] = {
212  "begin",
213  "parallel",
214  "connection",
215  "number" "of" "connections",
216  "element",
217  "node",
218  "force",
219  "rigid" "body",
220  "joint",
221  "beam",
222  "plate",
223  "rotor",
224  "aeromodal",
225  "aerodynamic" "element",
226  "electric" "bulk",
227  "electric",
228  "thermal",
229  "genel",
230  "hydraulic",
231  "bulk",
232  "loadable",
233  "driven",
234  "abstract",
235  "structural",
236  "parameter",
237  "weights",
238  "partition",
239  "partitioner",
240  "end",
241  NULL
242  };
243 
244  /* enum delle parole chiave */
245  enum KeyWords {
246  UNKNOWN = -1,
247  BEGIN = 0,
248  PARALLEL,
249  CONNECTION,
250  NUMBEROFCONNECTIONS,
251  ELEMENT,
252  NODE,
253  FORCE,
254  RIGIDBODY,
255  JOINT,
256  BEAM,
257  PLATE,
258  ROTOR,
259  AEROMODAL,
260  AERODYNAMICELEMENT,
261  ELECTRICBULK,
262  ELECTRIC,
263  THERMAL,
264  GENEL,
265  HYDRAULIC,
266  BULK,
267  LOADABLE,
268  DRIVEN,
269  ABSTRACT,
270  STRUCTURAL,
271  PARAMETER,
272  WEIGHTS,
273  PARTITION,
274  PARTITIONER,
275  END,
277  };
278 
279  /* tabella delle parole chiave */
280  KeyTable K(HP, sKeyWords);
281 
282  /* legge la distribuzione degli elementi sulle diverse CPU */
283 
284  try {
285  if (KeyWords(HP.GetDescription()) != BEGIN) {
286  pedantic_cerr("no explicit connections declared "
287  "for this input file" << std::endl);
288  return;
289  }
290 
291  } catch (EndOfFile) {
292  pedantic_cerr("no explicit connections declared "
293  "for this input file" << std::endl);
294  return;
295  }
296 
297  int iNumElems = 0;
298  int iNumNodes = 0;
299  if (KeyWords(HP.GetWord()) != PARALLEL) {
300  silent_cerr("Error: \"begin: parallel;\" expected at line "
301  << HP.GetLineData() << std::endl);
303  }
304 
305  while (true) {
306  switch (KeyWords(HP.GetDescription())) {
307  case WEIGHTS:
308  if (!HP.IsArg()) {
309  silent_cerr("Error: Weight flag expected "
310  "at line " << HP.GetLineData()
311  << std::endl);
313  }
314 
315  while (HP.IsArg()) {
316  if (HP.IsKeyWord("none")) {
318 
319  } else if (HP.IsKeyWord("vertices")
320  || HP.IsKeyWord("computation"))
321  {
323 
324  } else if (HP.IsKeyWord("communication")) {
325  wgtflag |= WEIGHT_COMM;
326 
327  } else if (HP.IsKeyWord("edges")) {
329 
330  } else {
331  silent_cerr("invalid weight "
332  " at line " << HP.GetLineData()
333  << std::endl);
335  }
336  }
337  break;
338 
339  case PARTITION:
341  for (int i = 0; i < iTotVertices; i++) {
342  if (!HP.IsArg()) {
343  silent_cerr("the partition "
344  "assignment is not complete; "
345  "only " << i << " vertices "
346  "input so far "
347  "at line " << HP.GetLineData()
348  << std::endl);
350  }
351  pParAmgProcs[i] = HP.GetInt();
352  if (pParAmgProcs[i] < 0 ||
354  {
355  silent_cerr("illegal value "
356  << pParAmgProcs[i]
357  << " for partition "
358  "assignment[" << i << "] "
359  "at line " << HP.GetLineData()
360  << "; must be between 0 and "
361  << DataCommSize - 1 << std::endl);
363  }
364  }
365  break;
366 
367  case PARTITIONER:
368  if (HP.IsKeyWord("metis")) {
369 #ifdef USE_METIS
371 #else /* ! USE_METIS */
372  silent_cerr("METIS partitioner not available; aborting..." << std::endl);
374 #endif /* ! USE_METIS */
375 
376  } else if (HP.IsKeyWord("chaco")) {
377 #ifdef USE_CHACO
379 #else /* ! USE_CHACO */
380  silent_cerr("CHACO partitioner not available; aborting..." << std::endl);
382 #endif /* ! USE_CHACO */
383  } else if (HP.IsKeyWord("manual")) {
385 
386  } else {
387  silent_cerr("unknown partitioner "
388  "at line " << HP.GetLineData()
389  << std::endl);
391  }
392  break;
393 
394  case END:
395  if (KeyWords(HP.GetWord()) != PARALLEL) {
396  silent_cerr("Error: \"end: parallel;\" expected "
397  "at line " << HP.GetLineData()
398  << "; aborting..." << std::endl);
400  }
401  goto endcycle;
402 
403  default:
404  silent_cerr("Unknown input at line "
405  << HP.GetLineData() << "; aborting..."
406  << std::endl);
408 
409  case NUMBEROFCONNECTIONS:
413 
414  Elem::Type CurrElType;
415  Node::Type CurrNdType;
416  unsigned int j = 0;
417 
418  for (int i = 0; i < iTotalExpConnections; i++) {
419  if (KeyWords(HP.GetDescription()) != CONNECTION) {
420  silent_cerr("Error: <Connection> "
421  "expected at line "
422  << HP.GetLineData()
423  << "; aborting..."
424  << std::endl);
426  }
427 
428  for (int k = 0; k < 2; k++) {
429  switch (KeyWords(HP.GetWord())) {
430  case ELEMENT:
431  switch (KeyWords(HP.GetWord())) {
432  case FORCE:
433  CurrElType = Elem::FORCE;
434  break;
435 
436  case RIGIDBODY:
437  CurrElType = Elem::BODY;
438  break;
439 
440  case JOINT:
441  CurrElType = Elem::JOINT;
442  break;
443 
444  case BEAM:
445  CurrElType = Elem::BEAM;
446  break;
447 
448  case PLATE:
449  CurrElType = Elem::PLATE;
450  break;
451 
452  case ROTOR:
453  CurrElType = Elem::INDUCEDVELOCITY;
454  break;
455 
456  case AEROMODAL:
457  CurrElType = Elem::AEROMODAL;
458  break;
459 
460  case AERODYNAMICELEMENT:
461  CurrElType = Elem::AERODYNAMIC;
462  break;
463 
464  case ELECTRICBULK:
465  CurrElType = Elem::ELECTRICBULK;
466  break;
467 
468  case ELECTRIC:
469  CurrElType = Elem::ELECTRIC;
470  break;
471 
472  case THERMAL:
473  CurrElType = Elem::THERMAL;
474  break;
475 
476  case GENEL:
477  CurrElType = Elem::GENEL;
478  break;
479 
480  case HYDRAULIC:
481  CurrElType = Elem::HYDRAULIC;
482  break;
483 
484  case BULK:
485  CurrElType = Elem::BULK;
486  break;
487 
488  case LOADABLE:
489  CurrElType = Elem::LOADABLE;
490  break;
491 
492  case DRIVEN:
493  CurrElType = Elem::DRIVEN;
494  break;
495 
496  default:
497  silent_cerr("Error: invalid element type "
498  "at line " << HP.GetLineData()
499  << "; aborting..." << std::endl);
501  }
502 
503  if (HP.IsArg()) {
504  j = HP.GetInt();
505 
506  } else {
507  silent_cerr("Error: label expected "
508  "for " << psElemNames[CurrElType]
509  << " element type at line "
510  << HP.GetLineData()
511  << "; aborting..." << std::endl);
513  }
514 
515  ppExpCntElems[iNumElems] =
516  *(ppFindElem(CurrElType, j));
517  if (ppExpCntElems[iNumElems] == NULL) {
518  silent_cerr("Error at line "
519  << HP.GetLineData() << ": "
520  << psElemNames[CurrElType]
521  << "(" << j << ") undefined; "
522  "aborting..." << std::endl);
524  }
525  iNumElems++;
526  break;
527 
528  case NODE:
529  switch (KeyWords(HP.GetWord())) {
530  case ABSTRACT:
531  CurrNdType = Node::ABSTRACT;
532  break;
533 
534  case STRUCTURAL:
535  CurrNdType = Node::STRUCTURAL;
536  break;
537 
538  case ELECTRIC:
539  CurrNdType = Node::ELECTRIC;
540  break;
541 
542  case THERMAL:
543  CurrNdType = Node::THERMAL;
544  break;
545 
546  case PARAMETER:
547  CurrNdType = Node::PARAMETER;
548  break;
549 
550  case HYDRAULIC:
551  CurrNdType = Node::HYDRAULIC;
552  break;
553 
554  default:
555  silent_cerr("Error: invalid node type "
556  "at line " << HP.GetLineData()
557  << "; aborting..." << std::endl);
559  }
560 
561  if (HP.IsArg()) {
562  j = HP.GetInt();
563 
564  } else {
565  silent_cerr("Error: label expected "
566  "at line " << HP.GetLineData()
567  << "; aborting..." << std::endl);
569  }
570 
571  ppExpCntNodes[iNumNodes] =
572  pFindNode(CurrNdType, j);
573  if (ppExpCntNodes[iNumNodes] == NULL ) {
574  silent_cerr("Error: at line "
575  << HP.GetLineData() << ":"
576  << psNodeNames[CurrNdType]
577  << "(" << j << ") undefined; "
578  << "aborting..." << std::endl);
580  }
581  iNumNodes++;
582  break;
583 
584  default:
585  silent_cerr("Unknown input at line "
586  << HP.GetLineData()
587  << "; aborting..." << std::endl);
589  }
590  }
591  }
592 
593  ASSERT(iNumNodes == iNumElems);
594  if (iNumNodes + iNumElems != 2*iTotalExpConnections) {
595  silent_cerr("Error: total number of nodes and elements"
596  " in the parallel section at line "
597  << HP.GetLineData()
598  << " is not consistent; aborting..." << std::endl);
600  }
601  break;
602  }
603  }
604 endcycle:;
605 }
606 /* Costruttore - End */
607 
608 /* Distruttore - begin */
610 {
611  DEBUGCOUTFNAME("SchurDataManager::~SchurDataManager");
612 
613  if (ppMyNodes != NULL) {
615  }
616 
617  if (ppMyElems != NULL) {
619  }
620 
621  if (ppMyIntElems != NULL) {
623  }
624 
625  if (ppIntNodes != NULL) {
627  }
628 
629  if (pLocalDofs != NULL) {
631  }
632 
633  if (pLocalIntDofs != NULL) {
635  }
636 
637  if (ppExpCntNodes != NULL) {
639  }
640 
641  if (ppExpCntElems != NULL) {
643  }
644 }
645 
646 integer
647 SchurDataManager::HowManyDofs(DofType who) const
648 {
649  switch(who) {
650  case TOTAL:
651  return iTotDofs;
652 
653  case LOCAL:
654  return iNumLocDofs;
655 
656  case INTERNAL:
657  return iNumIntDofs;
658 
659  case MYINTERNAL:
660  return iNumMyInt;
661 
662  default:
663  silent_cerr("SchurDataManager::HowManyDofs: "
664  "illegal request (" << unsigned(who) << ")"
665  << std::endl);
667  }
668 }
669 
670 integer *
671 SchurDataManager::GetDofsList(DofType who) const
672 {
673  switch(who) {
674  case LOCAL:
675  return pLocalDofs;
676 
677  case INTERNAL:
678  return pLocalIntDofs;
679 
680  case MYINTERNAL:
681  return pMyIntDofs;
682 
683  default:
684  silent_cerr("SchurDataManager::GetDofsList: "
685  "illegal request (" << unsigned(who) << ")"
686  << std::endl);
688  }
689 }
690 
691 Dof*
693 {
694  return pDofs;
695 }
696 
697 /* Ripartisce i nodi fra i processi e all'interno di ogni singolo processo */
698 void
700 {
701  DEBUGCOUT("Entering SchurDataManager::CreatePartition()" << std::endl);
702 
703  std::vector<std::vector<int> > adj(iTotVertices); /* adjacency */
704  std::vector<int> CommWgts(iTotVertices); /* communication weights */
705 
706  Adjacency Vertices; /* Struttura contenente le connessioni fra i vertici */
707  std::vector<int> VertexWgts(iTotVertices); /* Pesi dei vertici = dofs x ogni v. utile per METIS */
708  Vertices.pXadj = 0;
709  Vertices.pAdjncy = 0;
710  integer iMax = 0;
711  integer iRMax = 0;
712  int iCount = 0;
713  int iNumRt = 0;
714 
715  ASSERT(iTotVertices > 0);
716  ASSERT(DataCommSize > 0);
717 
718  /* Costruisco e inizializzo le due strutture */
719  SAFENEWARR(Vertices.pXadj, int, iTotVertices + 1);
720  memset(Vertices.pXadj, 0, (iTotVertices + 1)*sizeof(int));
721 
722  /* Ciclo per la scrittura degli array delle connessioni.
723  * Il ciclo viene ripetuto se per un vertice si ha un numero
724  * di connessioni superiore al max consentito per default
725  * iDefaultMaxConnectionsPerVertex */
726  int iMaxConnectionsPerVertex =
727  (iTotVertices < iDefaultMaxConnectionsPerVertex) ? iTotVertices
729  int GravityPos = 0, AirPropPos = 0;
730  int* pRotPos = NULL;
731  unsigned int* pRotLab = NULL;
732  if (!ElemData[Elem::INDUCEDVELOCITY].ElemContainer.empty()) {
734  SAFENEWARR(pRotPos, int, iNum);
735  SAFENEWARR(pRotLab, unsigned int, iNum);
736  memset(pRotPos, 0, iNum*sizeof(int));
737  memset(pRotLab, 0, iNum*sizeof(unsigned int));
738  }
739  SAFENEWARR(pLabelsList, unsigned int, iTotNodes);
740  memset(pLabelsList, 0, iTotNodes*sizeof(unsigned int));
741 
742  std::vector<const Node *> connectedNodes;
743 
744  while (true) {
745  InitList(Vertices.pXadj, iTotVertices + 1, 0);
746 
747  SAFENEWARR(Vertices.pAdjncy, int, iTotVertices*iMaxConnectionsPerVertex);
748  InitList(Vertices.pAdjncy, iTotVertices*iMaxConnectionsPerVertex, ADJ_UNDEFINED);
749  ASSERT(Elems.begin() != Elems.end());
750 
751  for (unsigned int i = 0; i < iTotNodes; i++) {
752  pLabelsList[i] = ppNodes[i]->GetLabel();
753  }
754 
755  /* ciclo sugli elementi per assemblare la struttura delle connessioni */
756  Node** ppCurrNode = NULL;
757  /* per numerare i nodi prendo la posizione del puntatore
758  * al nodo nell'array ppNodes */
759  int position;
760  iCount = iTotNodes;
761  iNumRt = 0;
762 
763  for (ElemVecType::const_iterator pTmpEl = Elems.begin();
764  pTmpEl != Elems.end();
765  pTmpEl++, iCount++)
766  {
767  if ((*pTmpEl)->GetElemType() == Elem::GRAVITY) {
768  GravityPos = iCount - iTotNodes;
769 
770  } else if ((*pTmpEl)->GetElemType() == Elem::AIRPROPERTIES) {
771  AirPropPos = iCount - iTotNodes;
772 
773  } else if ((*pTmpEl)->GetElemType() == Elem::INDUCEDVELOCITY) {
774  pRotPos[iNumRt] = iCount - iTotNodes;
775  pRotLab[iNumRt] = (*pTmpEl)->GetLabel();
776  iNumRt++;
777  }
778 
779  (*pTmpEl)->GetConnectedNodes(connectedNodes);
780 
781  /* peso dell'elemento */
782  integer dimA, dimB;
783  (*pTmpEl)->WorkSpaceDim(&dimA, &dimB);
784  VertexWgts[iCount] = dimA * dimB;
785 
786  CommWgts[iCount] = (*pTmpEl)->iGetNumDof()*(*pTmpEl)->iGetNumDof();
787 
788  for (std::vector<const Node *>::const_iterator i = connectedNodes.begin();
789  i != connectedNodes.end();
790  i++)
791  {
792  Vertices.pXadj[iCount + 1]++;
793 
794  /* trovo la pos. del nodo nella lista dei puntatori ai nodi */
795  Node::Type type = (*i)->GetNodeType();
796  unsigned label = (*i)->GetLabel();
797  ppCurrNode = SearchNode(NodeData[type].ppFirstNode,
798  NodeData[type].iNum, label);
799  position = ppCurrNode - ppNodes;
800 
801  /* Aggiungo al peso dell'elemento il numero di dofs
802  * di ciascun nodo connesso */
803 
804 #if 0 /* FIXME: i nodi hanno peso comp. nullo */
805  VertexWgts[iCount] += (*ppCurrNode)->iGetNumDof();
806 #endif
807 
808  /* aggiungo fra le connessioni dell'elemento il nodo attuale */
809  if ((iCount*iMaxConnectionsPerVertex) + Vertices.pXadj[iCount + 1] - 1 < iTotVertices*iMaxConnectionsPerVertex) {
810  Vertices.pAdjncy[(iCount*iMaxConnectionsPerVertex)
811  + Vertices.pXadj[iCount + 1] - 1] = position;
812  }
813 
814  /* aggiungo alle connessioni del nodo l'elemento attuale */
815  Vertices.pXadj[position + 1]++;
816  if ((position*iMaxConnectionsPerVertex) + Vertices.pXadj[position + 1] - 1 < iTotVertices*iMaxConnectionsPerVertex) {
817  Vertices.pAdjncy[(position*iMaxConnectionsPerVertex)
818  + Vertices.pXadj[position + 1] - 1] = iCount;
819  }
820 
821  /* peso (di comunicazione) del nodo */
822  CommWgts[position] = (*ppCurrNode)->iGetNumDof();
823  }
824  }
825 
826  if (iTotalExpConnections != 0) {
827  for (int i = 0; i < iTotalExpConnections; i++) {
828  int iNdPos, iElPos;
829  int j;
830 
831  for (j = 0; ppExpCntNodes[i] != ppNodes[j]; j++) {
832  NO_OP;
833  }
834  iNdPos = j;
835 
836  for (j = 0; ppExpCntElems[i] != Elems[j]; j++) {
837  NO_OP;
838  }
839  iElPos = j;
840 
841  Vertices.pXadj[iNdPos + 1]++;
842  Vertices.pXadj[iTotNodes + iElPos + 1]++;
843  Vertices.pAdjncy[(iNdPos*iMaxConnectionsPerVertex)
844  + Vertices.pXadj[iNdPos + 1] - 1] = iElPos + iTotNodes;
845  Vertices.pAdjncy[((iTotNodes + iElPos)*iMaxConnectionsPerVertex)
846  + Vertices.pXadj[iTotNodes + iElPos + 1] - 1] = iNdPos;
847  }
848  }
849 
850  iMax = 0;
851  for (int i = 0; i < iTotVertices; i++) {
852  if (Vertices.pXadj[i] > iMaxConnectionsPerVertex) {
853  iMax = Vertices.pXadj[i];
854  }
855  }
856 
857  if (iMax <= iMaxConnectionsPerVertex) {
858  break;
859  }
860 
861  iMaxConnectionsPerVertex = iMax;
862  SAFEDELETEARR(Vertices.pAdjncy);
863  Vertices.pAdjncy = 0;
864  }
865 
866  for (int i = 1; i <= iTotVertices; i++) {
867  Vertices.pXadj[i] += Vertices.pXadj[i - 1];
868  }
869 
870  /* Compatta il vettore delle adiacenze */
871  Pack(Vertices.pAdjncy, iTotVertices*iMaxConnectionsPerVertex);
872 
873  /* Chiamo la routine per ottere la partizione fra i diversi processi METIS.
874  * Se ne usano due diverse a seconda della dimensione della partizione */
875  if (pParAmgProcs == NULL) {
876  SAFENEWARR(pParAmgProcs, int, iTotVertices);
877 
878 #ifdef DEBUG
879  if (MyRank == 0) {
880  std::ofstream ofPartition;
881 
882  ofPartition.open("partition.debug");
883  ofPartition << "# METIS-like Input File" << std::endl
884  << "# Column 1: computational weight" << std::endl
885  << "# Column 2: communication weight" << std::endl
886  << "# Total Vertices: " << iTotVertices << std::endl
887  << "# Nodes" << std::endl;
888  for (unsigned int i = 0; i < iTotNodes; i++) {
889  ofPartition << "# " << i << " Node Type: "
890  << "(" << psNodeNames[ppNodes[i]->GetNodeType()] << ")"
891  << " Label: " << ppNodes[i]->GetLabel() << std::endl
892  << VertexWgts[i] << " " << CommWgts[i];
893  for (int j = Vertices.pXadj[i]; j < Vertices.pXadj[i + 1]; j++) {
894  ofPartition << " " << Vertices.pAdjncy[j];
895  }
896  ofPartition << std::endl;
897  }
898  ofPartition << "# Elements" << std::endl;
899  for (unsigned int i = 0; i < iTotElem; i++) {
900  ofPartition << "# " << i + iTotNodes << " Element Type: "
901  << "(" << psElemNames[Elems[i]->GetElemType()] << ")"
902  << " Label: " << Elems[i]->GetLabel() << std::endl
903  << VertexWgts[i + iTotNodes]
904  << " " << CommWgts[i + iTotNodes];
905  for (int j = Vertices.pXadj[i + iTotNodes];
906  j < Vertices.pXadj[i + iTotNodes + 1];
907  j++) {
908  ofPartition << " " << Vertices.pAdjncy[j];
909  }
910  ofPartition << std::endl;
911  }
912  ofPartition.close();
913  }
914 #endif /* DEBUG */
915 
916  switch (Partitioner) {
917  case PARTITIONER_MANUAL:
918  ASSERT(0);
919  break;
920 
921  case PARTITIONER_CHACO: {
922  int *vwgt = 0;
923  int *cwgt = 0;
924  int *ewgt = 0;
925 
926  if (wgtflag & WEIGHT_VERTICES) {
927  vwgt = &VertexWgts[0];
928  }
929 
930  if (wgtflag & WEIGHT_COMM) {
931  /* unsupported */
932  silent_cout("communication weights currently unsupported by CHACO; trying to emulate..." << std::endl);
933  cwgt = &CommWgts[0];
934  }
935 
936  if (wgtflag & WEIGHT_EDGES) {
937  /* unsupported ... */
938  silent_cout("edges weights currently unsupported by MBDyn" << std::endl);
939  }
940 
941  chaco_interface(iTotVertices,
942  Vertices.pXadj,
943  Vertices.pAdjncy,
944  vwgt,
945  cwgt,
946  ewgt,
947  DataCommSize,
948  pParAmgProcs);
949  break;
950  }
951 
952  case PARTITIONER_METIS: {
953  int *vwgt = 0;
954  int *cwgt = 0;
955 
956  if (wgtflag & WEIGHT_VERTICES) {
957  vwgt = &VertexWgts[0];
958  }
959 
960  if (wgtflag & WEIGHT_EDGES) {
961  /* unsupported ... */
962  silent_cout("edges weights currently unsupported by MBDyn" << std::endl);
963  }
964 
965  if (wgtflag & WEIGHT_COMM) {
966  cwgt = &CommWgts[0];
967  }
968 
969  mbdyn_METIS_PartGraph(iTotVertices,
970  Vertices.pXadj,
971  Vertices.pAdjncy,
972  vwgt,
973  cwgt,
974  NULL, /* ewgt */
975  DataCommSize,
976  pParAmgProcs);
977  break;
978  }
979 
980  default:
981  silent_cerr("no partition library is available; "
982  "partition assignments must be provided "
983  "manually" << std::endl);
985  }
986  }
987 
988  /* Stima del # vertici per blocco */
989  int MyDim = iTotVertices/DataCommSize + DataCommSize;
990 
991  /* Lista dei nodi appartenenti a questo processo */
992  Adjacency InterfNodes; /* nodi di interfaccia */
993 
994  /* Anche qui si usa un ciclo while per verificare
995  * se il numero di nodi di interfaccia per processo
996  * ipotizzato, iMaxInterfNodes, sufficiente */
997  integer iMaxInterfNodes = iDefaultMaxInterfNodes;
998  InterfNodes.pXadj = NULL;
999  InterfNodes.pAdjncy = NULL;
1000  SAFENEWARR(InterfNodes.pXadj, int, DataCommSize);
1001  SAFENEWARR(ppMyNodes, Node*, 2*MyDim);
1002  memset(InterfNodes.pXadj, 0, DataCommSize*sizeof(int));
1003  memset(ppMyNodes, 0, 2*MyDim*sizeof(Node*));
1004 
1005  while (true) {
1006  InitList(InterfNodes.pXadj, DataCommSize, 0);
1007  SAFENEWARR(InterfNodes.pAdjncy, int, DataCommSize*iMaxInterfNodes*8);
1008  memset(InterfNodes.pAdjncy, 0, DataCommSize*iMaxInterfNodes*8*sizeof(int));
1009  InitList(InterfNodes.pAdjncy, DataCommSize*iMaxInterfNodes*8, ADJ_UNDEFINED);
1010 
1011  iNumLocNodes = 0;
1012  iNumLocDofs = 0;
1013  iCount = 0;
1014  int TmpPrc = 0;
1015  bool bIsNInterf;
1016 
1017  for (unsigned int i = 0; i < iTotNodes; i++) {
1018  bIsNInterf = true;
1019 
1020  if (pParAmgProcs[i] == MyRank) {
1021  /* se uno dei nodi e' connesso ad un elemento non appartenente
1022  * a questo processo e' un nodo di interfaccia */
1023  for (int j = Vertices.pXadj[i]; j < Vertices.pXadj[i + 1]; j++) {
1024  TmpPrc = pParAmgProcs[Vertices.pAdjncy[j]];
1025  if (TmpPrc != MyRank) {
1026  InterfNodes.pXadj[TmpPrc] += 1;
1027  int k = TmpPrc * iMaxInterfNodes * 2 + InterfNodes.pXadj[TmpPrc] - 1;
1028  if (k < DataCommSize * iMaxInterfNodes * 2) {
1029  InterfNodes.pAdjncy[k] = i;
1030  bIsNInterf = false;
1031  }
1032  }
1033  }
1034 
1035  if (bIsNInterf) {
1036  ppMyNodes[iCount] = ppNodes[i];
1037  iNumLocNodes++;
1038  iNumLocDofs += ppMyNodes[iCount]->iGetNumDof();
1039  iCount++;
1040  }
1041  }
1042  }
1043 
1044  iMax = 0;
1045  for (int i = 0; i < DataCommSize; i++) {
1046  if (InterfNodes.pXadj[i] > iMaxInterfNodes) {
1047  iMax = InterfNodes.pXadj[i];
1048  }
1049  }
1050 
1051  DataComm.Allreduce(&iMax, &iRMax, 1, MPI::INT, MPI::MAX);
1052  iMax = iRMax;
1053  if (iMax <= iMaxInterfNodes) {
1054  break;
1055  }
1056 
1057  iMaxInterfNodes = iMax;
1058  SAFEDELETEARR(InterfNodes.pAdjncy);
1059  InterfNodes.pAdjncy = 0;
1060  }
1061 
1062  /* Scambio i dati riguardo ai nodi di interfaccia */
1063  MPI::Request* pRReq = NULL;
1064  MPI::Request* pSReq = NULL;
1065 
1066  SAFENEWARRNOFILL(pRReq, MPI::Request, DataCommSize);
1067  SAFENEWARRNOFILL(pSReq, MPI::Request, DataCommSize);
1068  for (int i = 0; i < DataCommSize; i++) {
1069  pRReq[i] = MPI::REQUEST_NULL;
1070  pSReq[i] = MPI::REQUEST_NULL;
1071  }
1072 
1073  const int DIM_TAG = 10;
1074 
1075  for (int i = 0; i < DataCommSize; i++) {
1076 #if 0
1077  if (i != MyRank)
1078 #endif
1079  {
1080  pRReq[i] = DataComm.Irecv(InterfNodes.pAdjncy + iMaxInterfNodes + i*iMaxInterfNodes*2,
1081  iMaxInterfNodes, MPI::INT, i, DIM_TAG);
1082  pSReq[i] = DataComm.Isend(InterfNodes.pAdjncy + i*iMaxInterfNodes*2,
1083  iMaxInterfNodes, MPI::INT, i, DIM_TAG);
1084  }
1085  }
1086 
1087  /* lista degli elementi appartenenti a questo processo */
1088  SAFENEWARR(ppMyElems, Elem*, 2*MyDim);
1089  memset(ppMyElems, 0, 2*MyDim*sizeof(Elem*));
1090 
1091  /* Trattamento elementi particolari */
1092  int move = 0;
1093 
1094  /* Gravity */
1095  if (!ElemData[Elem::GRAVITY].ElemContainer.empty()) {
1096  /* FIXME: there's a better way to find GravityPos and so on... */
1097  ppMyElems[iNumLocElems] = Elems[GravityPos];
1098  iNumLocElems++;
1099  pParAmgProcs[GravityPos + iTotNodes] = -1;
1100  move++;
1101  }
1102 
1103  /* Air Properties */
1104  if (!ElemData[Elem::AIRPROPERTIES].ElemContainer.empty()) {
1105  ppMyElems[iNumLocElems] = Elems[AirPropPos];
1106  iNumLocElems++;
1107  pParAmgProcs[AirPropPos + iTotNodes] = -1;
1108  move++;
1109  }
1110 
1111  /* Induced Velocity elements */
1112  int iMyTotRot = 0;
1113  integer* pMyRot = NULL;
1114  integer iIVIsMine = 0;
1115  if (iNumRt != 0) {
1116  SAFENEWARR(pMyRot, integer, iNumRt);
1117  for (ElemMapType::const_iterator i = ElemData[Elem::AERODYNAMIC].ElemContainer.begin();
1119  i++)
1120  {
1121  const AerodynamicElem *pAero = dynamic_cast<AerodynamicElem *>(i->second);
1122  ASSERT(pAero != NULL);
1123  const InducedVelocity *pIV = pAero->pGetInducedVelocity();
1124 
1125  if (pIV != NULL) {
1126  int pos = 0;
1127  unsigned int pTmpLab = pIV->GetLabel();
1128 
1129  for (int k = 0; k < iNumRt; k++) {
1130  if (pTmpLab == pRotLab[k]) {
1131  pos = pRotPos[k];
1132  }
1133  }
1134 
1135  for (int j = 0; j < iMyTotRot; j++) {
1136  if (pos != pMyRot[j]) {
1137  pMyRot[iMyTotRot] = pos;
1138  iMyTotRot++;
1139  }
1140  }
1141 
1142  if (!iMyTotRot) {
1143  pMyRot[iMyTotRot] = pos;
1144  iMyTotRot++;
1145  }
1146  }
1147  }
1148 
1149  /* Costruisco i communicators per i rotori */
1150  SAFENEWARRNOFILL(pIndVelComm, MPI::Intracomm, iNumRt);
1151 
1152  int color, key;
1153  for (int i = 0; i < iNumRt; i++) {
1154  if (iMyTotRot == 0) {
1155  color = MPI::UNDEFINED;
1156  key = MyRank;
1157  } else {
1158  color = 1;
1159  if (pParAmgProcs[pMyRot[i] + iTotNodes] == MyRank) {
1160  silent_cout("InducedVelocity(" << Elems[pRotPos[i]]->GetLabel() << ")"
1161  << " assigned to process "
1162  << MyRank << std::endl);
1163  iIVIsMine = 1;
1164  key = 0;
1165  } else {
1166  key = MyRank + 1;
1167  }
1168  }
1169  pIndVelComm[i] = MBDynComm.Split(color, key);
1170 #if 0
1171  IndVelComm[i] = MPI::COMM_WORLD.Split(color, key);
1172 #endif
1173  InducedVelocity *r = dynamic_cast<InducedVelocity *>(Elems[pRotPos[i]]);
1174  ASSERT(r != 0);
1175  r->InitializeIndVelComm(pIndVelComm + i);
1176  }
1177 
1178  for (int i = 0; i < iMyTotRot; i++) {
1179  ppMyElems[iNumLocElems] = Elems[pMyRot[i]];
1180  if (pParAmgProcs[pMyRot[i] + iTotNodes] == MyRank) {
1182  } else {
1183  move++;
1184  }
1185  iNumLocElems++;
1186  pParAmgProcs[pMyRot[i] + iTotNodes] = -1;
1187  }
1188  }
1189 
1190  for (unsigned int i = 0; i < iTotVertices - iTotNodes; i++) {
1191  if (pParAmgProcs[iTotNodes + i] == MyRank) {
1194  iNumLocElems++;
1195  }
1196  }
1197 
1198  /* initialize local element iterator */
1200 
1201 #ifdef USE_MPI_WAIT
1202  MPI::Request::Waitall(DataCommSize, pRReq);
1203  MPI::Request::Waitall(DataCommSize, pSReq);
1204 #else /* ! USE_MPI_WAIT */
1205  /* Verifico la ricezione dei nodi di interfaccia */
1206  bool bRecvFlag = false, bSentFlag = false;
1207  while (true) {
1208  if (!bRecvFlag) {
1209  bRecvFlag = MPI::Request::Testall(DataCommSize, pRReq);
1210  }
1211 
1212  if (!bSentFlag) {
1213  bSentFlag = MPI::Request::Testall(DataCommSize, pSReq);
1214  }
1215 
1216  if (bRecvFlag && bSentFlag) {
1217  break;
1218  }
1219 
1220  MYSLEEP(1000);
1221  }
1222 #endif /* ! USE_MPI_WAIT */
1223 
1224  /* ordino i nodi interfaccia */
1225  std::sort(InterfNodes.pAdjncy,
1226  InterfNodes.pAdjncy + iMaxInterfNodes * 2 * DataCommSize);
1227 
1228  /* elimino le ripetizioni */
1229  int* p = std::unique(InterfNodes.pAdjncy,
1230  InterfNodes.pAdjncy + iMaxInterfNodes * 2 * DataCommSize);
1231 
1232  /* dimensione effettiva dell'interfaccia locale
1233  * il -1 e' dovuto al primo valore che sempre pari a -1 */
1234  /*
1235  * InterfNodes.pAdjncy[0] == -1
1236  * InterfNodes.pAdjncy[iNumIntNodes] == -1
1237  */
1238  iNumIntNodes = p - &InterfNodes.pAdjncy[1];
1239 
1240  unsigned int* llabels = NULL;
1241  Node::Type* lTypes = NULL;
1242  SAFENEWARR(llabels, unsigned int, iNumIntNodes);
1245  for (int i = 0; i < iNumIntNodes; i++) {
1246  ppIntNodes[i] = ppNodes[InterfNodes.pAdjncy[i + 1]];
1247  llabels[i] = ppIntNodes[i]->GetLabel();
1248  lTypes[i] = ppIntNodes[i]->GetNodeType();
1249  }
1250 
1251  /* scrittura dei dofs degli stati interni di elementi che connettono
1252  * un nodo interno ad uno di interfaccia, sulla lista dei dof
1253  * di interfaccia */
1254  int* pPosIntElems = NULL;
1255  if (iNumLocElems != 0) {
1256  SAFENEWARR(pPosIntElems, int, iNumLocElems);
1258 
1259  memset(pPosIntElems, 0, iNumLocElems*sizeof(int));
1260  memset(ppMyIntElems, 0, iNumLocElems*sizeof(Elem*));
1261  }
1262 
1263  for (int i = 0; i < iNumLocElems; i++) {
1264  if (ppMyElems[i]->iGetNumDof() != 0) {
1265  Elem::Type CType = ppMyElems[i]->GetElemType();
1266  switch (CType) {
1267  case Elem::INDUCEDVELOCITY:
1268  if (iIVIsMine == 1) {
1269  ppMyElems[i]->GetConnectedNodes(connectedNodes);
1270  for (std::vector<const Node *>::const_iterator j = connectedNodes.begin();
1271  j != connectedNodes.end();
1272  j++)
1273  {
1274  unsigned int* p = std::find(llabels, llabels + iNumIntNodes, (*j)->GetLabel());
1275  if (p != llabels + iNumIntNodes) {
1277  iNumIntDofs += ppMyElems[i]->iGetNumDof();
1278  iNumLocDofs -= ppMyElems[i]->iGetNumDof();
1279  pPosIntElems[iNumIntElems] = i;
1280  iNumIntElems++;
1281  break;
1282  }
1283  }
1284  }
1285  break;
1286 
1287  case Elem::HYDRAULIC:
1288  case Elem::GENEL:
1289  case Elem::ELECTRIC:
1290  case Elem::THERMAL:
1291  break;
1292 
1293  default:
1294  ppMyElems[i]->GetConnectedNodes(connectedNodes);
1295  for (std::vector<const Node *>::const_iterator j = connectedNodes.begin();
1296  j != connectedNodes.end();
1297  j++)
1298  {
1299  unsigned int* p = std::find(llabels, llabels + iNumIntNodes, (*j)->GetLabel());
1300  if (p != llabels + iNumIntNodes) {
1302  iNumIntDofs += ppMyElems[i]->iGetNumDof();
1303  iNumLocDofs -= ppMyElems[i]->iGetNumDof();
1304  pPosIntElems[iNumIntElems] = i;
1305  iNumIntElems++;
1306  break;
1307  }
1308  }
1309  break;
1310  }
1311  }
1312  }
1313 
1314  /* determina la liste dei dofs locali ed adiacenti suddivisa per processi,
1315  * secondo la struttura Adjacency */
1317 
1318  iCount = 0;
1319  for (int i = 0; i < iNumLocNodes; i++) {
1320  if (ppMyNodes[i]->iGetNumDof() != 0) {
1321  integer First = (ppMyNodes[i])->iGetFirstIndex();
1322 
1323  pLocalDofs[iCount] = First + 1;
1324  iCount++;
1325  for (unsigned int j = 1; j < ppMyNodes[i]->iGetNumDof(); j++) {
1326  pLocalDofs[iCount] = First + j + 1;
1327  iCount++;
1328  }
1329  }
1330  }
1331 
1332  /* Aggiungo i dofs degli elementi locali */
1333  integer TmpDofNum;
1334  int i2Count = 0;
1335 
1336  for (int i = move; i < iNumLocElems; i++) {
1337  TmpDofNum = ppMyElems[i]->iGetNumDof();
1338  if (TmpDofNum != 0) {
1339  if (i != pPosIntElems[i2Count]) {
1340  ElemWithDofs* pWithDofs = dynamic_cast<ElemWithDofs *>(ppMyElems[i]);
1341  integer First = (pWithDofs)->iGetFirstIndex();
1342 
1343  pLocalDofs[iCount] = First + 1;
1344  iCount++;
1345  for (int j = 1; j < TmpDofNum; j++) {
1346  pLocalDofs[iCount] = First + 1 + j;
1347  iCount++;
1348  }
1349 
1350  } else {
1351  i2Count++;
1352  }
1353  }
1354  }
1355 
1356  /* scrivo ora la lista dei dofs interfaccia */
1357  for (int i = 1; i < iNumIntNodes + 1; i++) {
1358  iNumIntDofs += ppNodes[InterfNodes.pAdjncy[i]]->iGetNumDof();
1359  }
1360 
1363 
1364  iCount = 0;
1365  i2Count = 0;
1366  for (int i = 1; i <= iNumIntNodes; i++) {
1367  if (ppNodes[InterfNodes.pAdjncy[i]]->iGetNumDof() != 0) {
1368  integer First = ppNodes[InterfNodes.pAdjncy[i]]->iGetFirstIndex();
1369 
1370  pLocalIntDofs[iCount] = First + 1;
1371  iCount++;
1372  if (pParAmgProcs[InterfNodes.pAdjncy[i]] == MyRank) {
1373  pMyIntDofs[i2Count] = First + 1;
1374  i2Count++;
1375  }
1376 
1377  for (unsigned int j = 1;
1378  j < ppNodes[InterfNodes.pAdjncy[i]]->iGetNumDof();
1379  j++)
1380  {
1381  /* il - serve a distinguere questi dofs da quelli interni */
1382  pLocalIntDofs[iCount] = (First + 1 + j);
1383  iCount++;
1384  if (pParAmgProcs[InterfNodes.pAdjncy[i]] == MyRank) {
1385  pMyIntDofs[i2Count] = (First + 1 + j);
1386  i2Count++;
1387  }
1388  }
1389  }
1390  }
1391 
1392  /* Interfaccia degli elementi locali */
1393  for (int i = 0; i < iNumIntElems; i++) {
1394  TmpDofNum = ppMyIntElems[i]->iGetNumDof();
1395  ElemWithDofs* pWithDofs = dynamic_cast<ElemWithDofs *>(ppMyIntElems[i]);
1396  integer First = (pWithDofs)->iGetFirstIndex();
1397  pLocalIntDofs[iCount] = First + 1;
1398  iCount++;
1399  for (int j = 1; j < TmpDofNum; j++) {
1400  pLocalIntDofs[iCount] = First + 1 + j;
1401  iCount++;
1402  }
1403  }
1404 
1405  iNumMyInt = i2Count;
1406 
1407  if ( Vertices.pXadj != NULL) {
1408  SAFEDELETEARR(Vertices.pXadj);
1409  }
1410 
1411  if ( Vertices.pAdjncy != NULL) {
1412  SAFEDELETEARR(Vertices.pAdjncy);
1413  }
1414 
1415  if ( InterfNodes.pXadj != NULL) {
1416  SAFEDELETEARR(InterfNodes.pXadj);
1417  }
1418 
1419  if ( InterfNodes.pAdjncy != NULL) {
1420  SAFEDELETEARR( InterfNodes.pAdjncy);
1421  }
1422 
1423  if ( pParAmgProcs != NULL) {
1425  }
1426 
1427  if ( pLabelsList != NULL) {
1429  }
1430 
1431  if ( pPosIntElems != NULL) {
1432  SAFEDELETEARR(pPosIntElems);
1433  }
1434 
1435  if ( lTypes != NULL) {
1436  SAFEDELETEARR(lTypes);
1437  }
1438 
1439  if ( llabels != NULL) {
1440  SAFEDELETEARR(llabels);
1441  }
1442 
1443  OutputPartition();
1444 }
1445 
1446 void
1448 {
1449  silent_cout("Making partition file...");
1450 
1451  /* Inizializzazione */
1453 
1454  time_t tCurrTime(time(NULL));
1455  OutHdl.Partition()
1456  << "# Partition file for MBDyn. Time: " << ctime(&tCurrTime) << std::endl
1457  << "# Partition produced ";
1458  switch (Partitioner) {
1459  case PARTITIONER_MANUAL:
1460  OutHdl.Partition()
1461  << "manually";
1462  break;
1463 
1464  case PARTITIONER_METIS:
1465  OutHdl.Partition()
1466  << "with METIS";
1467  break;
1468 
1469  case PARTITIONER_CHACO:
1470  OutHdl.Partition()
1471  << "with CHACO";
1472  break;
1473 
1474  default:
1475  ASSERT(0);
1476  }
1477  OutHdl.Partition()
1478  << std::endl;
1479 
1480  /* Dati */
1481  OutHdl.Partition()
1482  << std::endl
1483  << std::endl
1484  << "# Control data to verify the partition "
1485  << std::endl
1486  << std::endl
1487  << "Total number of processes: " << DataCommSize << ";" << std::endl
1488  << "Process #: " << MyRank << ";" << std::endl
1489  << "Total number of Nodes: " << iTotNodes << ";" << std::endl
1490  << "Total number of Elements: " << Elems.size() << ";" << std::endl
1491  << "Total number of Dofs: " << iTotDofs << ";" << std::endl;
1492 
1493  OutHdl.Partition()
1494  << std::endl
1495  << std::endl
1496  << "Local Dofs number: " << iNumLocDofs << std::endl;
1497 
1498  OutHdl.Partition()
1499  << std::endl
1500  << "Local Nodes: " << iNumLocNodes << std::endl;
1501  for (int i = 0; i < iNumLocNodes; i++) {
1502  ASSERT(ppMyNodes[i] != NULL);
1503  OutHdl.Partition()
1504  << "Node Type: "
1505  << "(" << psNodeNames[(ppMyNodes[i])->GetNodeType()] << ")"
1506  << " Label: " << ppMyNodes[i]->GetLabel()
1507  << " Dofs #: " << ppMyNodes[i]->iGetNumDof()
1508  << std::endl;
1509  }
1510 
1511  OutHdl.Partition()
1512  << std::endl
1513  << std::endl
1514  << "Local Elements: "<< iNumLocElems << std::endl;
1515  for (int i = 0; i < iNumLocElems; i++) {
1516  ASSERT(ppMyElems[i] != NULL);
1517  OutHdl.Partition()
1518  << "Element Type: "
1519  << "(" << psElemNames[ppMyElems[i]->GetElemType()] << ")"
1520  << " Label: " << ppMyElems[i]->GetLabel()
1521  << " Dofs #: " << ppMyElems[i]->iGetNumDof()
1522  << std::endl;
1523  }
1524 
1525  OutHdl.Partition()
1526  << std::endl
1527  << std::endl
1528  << "Interface Dofs number: " << iNumIntDofs << std::endl
1529  << std::endl
1530  << "Interface Nodes: " << iNumIntNodes << std::endl;
1531  for (int i = 0; i < iNumIntNodes; i++) {
1532  ASSERT(ppIntNodes[i] != NULL);
1533  OutHdl.Partition()
1534  << "Node Type: "
1535  << "(" << psNodeNames[(ppIntNodes[i])->GetNodeType()] << ")"
1536  << " Label: " << ppIntNodes[i]->GetLabel()
1537  << " Dofs #: " << ppIntNodes[i]->iGetNumDof()
1538  << std::endl;
1539  }
1540 
1541  OutHdl.Partition()
1542  << std::endl
1543  << std::endl
1544  << "Elements whose internal dofs are interface dofs: "
1545  << iNumIntElems << std::endl;
1546  for (int i = 0; i < iNumIntElems; i++) {
1547  OutHdl.Partition()
1548  << "Element Type: "
1549  << "(" << psElemNames[(ppMyIntElems[i])->GetElemType()] << ")"
1550  << " Label: " << ppMyIntElems[i]->GetLabel()
1551  << " Dofs #: " << ppMyIntElems[i]->iGetNumDof()
1552  << std::endl;
1553  }
1554 
1555 #ifdef DEBUG
1556  OutHdl.Partition()
1557  << std::endl
1558  << std::endl
1559  << "Local Dofs List:" << std::endl;
1560 
1561  for (int i = 0; i < iNumLocDofs; i++) {
1562  OutHdl.Partition() << pLocalDofs[i] << " ";
1563  if (!(iNumLocDofs%10)) {
1564  OutHdl.Partition() << std::endl;
1565  }
1566  }
1567  if (iNumLocDofs%10) {
1568  OutHdl.Partition()
1569  << std::endl;
1570  }
1571 
1572  OutHdl.Partition()
1573  << std::endl
1574  << std::endl
1575  << "Local Interface Dofs List:" <<std::endl;
1576  for (int i = 0; i < iNumIntDofs; i++) {
1577  OutHdl.Partition() << pLocalIntDofs[i] << " ";
1578  if (!(i%10)) {
1579  OutHdl.Partition() << std::endl;
1580  }
1581  }
1582  OutHdl.Partition() << std::endl;
1583 #endif /* DEBUG */
1584 
1585  silent_cout("... partition file done" << std::endl);
1586 }
1587 
1588 /* compatta il vettore delle adiacenze */
1589 void
1590 SchurDataManager::Pack(int* pList, int dim)
1591 {
1592  int* pOld = pList;
1593  int* pNew = pList;
1594 
1595  for (; pOld < pList + dim; pOld++) {
1596  if (*pOld != ADJ_UNDEFINED) {
1597  *pNew = *pOld;
1598  pNew++;
1599  }
1600  }
1601 }
1602 
1603 /* Inizializza le varie liste di interi usate in Create Partition */
1604 void
1605 SchurDataManager::InitList(int* list, int dim, int value)
1606 {
1607  for (int i = 0; i < dim; i++) {
1608  list[i] = value;
1609  }
1610 }
1611 
1612 void
1613 SchurDataManager::InitList(float* list, int dim, int value)
1614 {
1615  for (int i = 0; i < dim; i++) {
1616  list[i] = value;
1617  }
1618 }
1619 
1620 
1621 /* Trova la posizione del nodo nell'array ppNodes */
1622 Node**
1623 SchurDataManager::SearchNode(Node** ppFirst, int dim, unsigned int& label)
1624 {
1625  unsigned int* pbegin = pLabelsList + (ppFirst - ppNodes);
1626  unsigned int* p = std::find(pbegin, pbegin + dim, label);
1627 
1628  return ppNodes + (p - pLabelsList);
1629 }
1630 
1631 
1632 void
1633 SchurDataManager::AssRes(VectorHandler& ResHdl, doublereal dCoef) throw(ChangedEquationStructure)
1634 {
1635  DEBUGCOUT("Entering SchurDataManager::AssRes()" << std::endl);
1636  ASSERT(pWorkVec != NULL);
1637 
1638  try {
1639  DataManager::AssRes(ResHdl, dCoef, MyElemIter, *pWorkVec);
1640 
1641  } catch (ChangedEquationStructure) {
1642  Elem *pEl = NULL;
1643 
1644  if (MyElemIter.bGetCurr(pEl) == true) {
1645  silent_cerr("Jacobian reassembly requested by "
1646  << psElemNames[pEl->GetElemType()]
1647  << "(" << pEl->GetLabel() << "); "
1648  "currently unsupported by Schur data manager."
1649  << std::endl);
1650 
1651  } else {
1652  silent_cerr("Jacobian reassembly requested "
1653  "by an element; currently unsupported "
1654  "by Schur data manager." << std::endl);
1655  }
1656 
1658  }
1659 }
1660 
1661 void
1663 {
1664  DEBUGCOUT("Entering SchurDataManager::AssJac()" << std::endl);
1665  ASSERT(pWorkMat != NULL);
1666 
1667  DataManager::AssJac(JacHdl, dCoef, MyElemIter, *pWorkMat);
1668 }
1669 /* End of AssJac */
1670 
1671 void
1672 SchurDataManager::Update(void) const
1673 {
1674  DEBUGCOUT("Entering SchurDataManager::Update()" << std::endl);
1675 
1676  /* Nodi */
1677  for (int i = 0; i < iNumLocNodes; i++) {
1678  ASSERT(ppMyNodes[i] != NULL);
1680  }
1681 
1682  /* Nodi di interfaccia */
1683  for (int i = 0; i < iNumIntNodes; i++) {
1684  ASSERT(ppIntNodes[i] != NULL);
1686  }
1687 
1688  /* Elementi */
1689  for (int i = 0; i < iNumLocElems; i++) {
1690  ASSERT(ppMyElems[i] != NULL);
1692  }
1693 }
1694 /* End of Update */
1695 
1696 void
1698 {
1699  DEBUGCOUT("Entering SchurDataManager::DerivativesUpdate()" << std::endl);
1700 
1701  /* Nodi */
1702  for (int i = 0; i < iNumLocNodes; i++) {
1703  ASSERT(ppMyNodes[i] != NULL);
1704  if (ppMyNodes[i]->GetNodeType() == Node::STRUCTURAL) {
1705  ((StructNode*)ppMyNodes[i])->DerivativesUpdate(*pXCurr, *pXPrimeCurr);
1706  } else {
1708  }
1709  }
1710 
1711  /* Nodi adiacenti i cui valori influenzano gli assemblaggi degli elementi */
1712  for (int i = 0; i < iNumIntNodes; i++) {
1713  ASSERT(ppIntNodes[i] != NULL);
1714  if (ppIntNodes[i]->GetNodeType() == Node::STRUCTURAL) {
1715  ((StructNode*)ppIntNodes[i])->DerivativesUpdate(*pXCurr, *pXPrimeCurr);
1716  } else {
1718  }
1719  }
1720 
1721  /* Elementi */
1722  for (int i = 0; i < iNumLocElems; i++) {
1723  ASSERT(ppMyElems[i] != NULL);
1725  }
1726 }
1727 /* End of DerivativeUpdate */
1728 
1729 
1730 void
1732  VectorHandler& XPrev, VectorHandler& XPPrev) const
1733 {
1734  DEBUGCOUT("Entering SchurDataManager::BeforePredict()" << std::endl);
1735 
1736  /* Nodi */
1737  for (int i = 0; i < iNumLocNodes; i++) {
1738  ASSERT(ppMyNodes[i] != NULL);
1739  ppMyNodes[i]->BeforePredict(X, XP, XPrev, XPPrev);
1740  }
1741 
1742  /* Nodi adiacenti i cui valori influenzano gli assemblaggi degli elementi */
1743  for (int i = 0; i < iNumIntNodes; i++) {
1744  ASSERT(ppIntNodes[i] != NULL);
1745  ppIntNodes[i]->BeforePredict(X, XP, XPrev, XPPrev);
1746  }
1747 
1748  /* Elementi */
1749  for (int i = 0; i < iNumLocElems; i++) {
1750  ASSERT(ppMyElems[i] != NULL);
1751  ppMyElems[i]->BeforePredict(X, XP, XPrev, XPPrev);
1752  }
1753 }
1754 /* End of BeforePredict */
1755 
1756 void
1758 {
1759  DEBUGCOUT("Entering SchurDataManager::AfterPredict()" << std::endl);
1760 
1761  /* Nodi */
1762  for (int i = 0; i < iNumLocNodes; i++) {
1763  ASSERT(ppMyNodes[i] != NULL);
1766  }
1767 
1768  /* Nodi adiacenti i cui valori influenzano gli assemblaggi degli elementi */
1769  for (int i = 0; i < iNumIntNodes; i++) {
1770  ASSERT(ppIntNodes[i] != NULL);
1773  }
1774 
1775  /* Elementi */
1776  for (int i = 0; i < iNumLocElems; i++) {
1777  ASSERT(ppMyElems[i] != NULL);
1780  }
1781 }
1782 /* End of AfterPredict */
1783 
1784 
1785 void
1787 {
1788  DEBUGCOUT("Entering SchurDataManager::AfterConvergence()" << std::endl);
1789 
1790  /* Nodi */
1791  for (int i = 0; i < iNumLocNodes; i++) {
1792  ASSERT(ppMyNodes[i] != NULL);
1795  }
1796 
1797  /* Nodi adiacenti i cui valori influenzano gli assemblaggi degli elementi */
1798  for (int i = 0; i < iNumIntNodes; i++) {
1799  ASSERT(ppIntNodes[i] != NULL);
1802  }
1803 
1804  /* Elementi */
1805  for (int i = 0; i < iNumLocElems; i++) {
1806  ASSERT(ppMyElems[i] != NULL);
1809  }
1810 
1811  /* Restart condizionato */
1812  switch (RestartEvery) {
1813  case ITERATIONS:
1815  iCurrRestartIter = 0;
1816  ((SchurDataManager*)this)->MakeRestart();
1817  }
1818  break;
1819 
1820  case TIME: {
1821  doublereal dT = dGetTime();
1822  if (dT - dLastRestartTime >= dRestartTime) {
1823  dLastRestartTime = dT;
1824  const_cast<SchurDataManager *>(this)->MakeRestart();
1825  }
1826  break;
1827  }
1828 
1829  default:
1830  break;
1831  }
1832 }
1833 /* End of AfterConvergence */
1834 
1835 
1836 /* stampa i risultati */
1837 void
1838 SchurDataManager::Output(bool force) const
1839 {
1840  DEBUGCOUT("Entering SchurDataManager::Output()" << std::endl);
1841 
1842  /* output only at multiples of iOutputFrequency */
1843  if ((!force) && !pOutputMeter->dGet()) {
1844  return;
1845  }
1846 
1847  /* Nodi */
1848  for (int i = 0; i < iNumLocNodes; i++) {
1849  ASSERT(ppMyNodes[i] != NULL);
1850  ppMyNodes[i]->Output(OutHdl);
1851  }
1852 
1853  /* Nodi di interfaccia */
1854  for (int i = 0; i < iNumIntNodes; i++) {
1855  ASSERT(ppIntNodes[i] != NULL);
1856  ppIntNodes[i]->Output(OutHdl);
1857  }
1858 
1859  /* Elementi */
1860  for (int i = 0; i < iNumLocElems; i++) {
1861  ASSERT(ppMyElems[i] != NULL);
1862  ppMyElems[i]->Output(OutHdl);
1863  }
1864 }
1865 
1866 /* End Output */
1867 
1868 /* End SchurDataManager */
1869 
1870 #endif /* USE_MPI */
1871 
eRestart RestartEvery
Definition: dataman.h:164
integer * GetDofsList(DofType who) const
integer * pLocalIntDofs
Definition: schurdataman.h:78
ElemVecType Elems
Definition: dataman.h:609
enum SchurDataManager::PartitionLibrary Partitioner
unsigned wgtflag
Definition: schurdataman.h:91
ElemContainerType ElemContainer
Definition: dataman.h:591
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Definition: node.h:67
#define DEBUGCOUTFNAME(fname)
Definition: myassert.h:256
void DerivativesUpdate(void) const
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
void Update(void) const
doublereal dGetTime(void) const
Definition: dataman2.cc:165
const integer iDefaultMaxInterfNodes
Definition: schurdataman.cc:66
OutputHandler OutHdl
Definition: dataman.h:105
unsigned int iTotNodes
Definition: dataman.h:748
static int position(const MathParser::MathArgs &args)
Definition: modelns.cc:99
Elem ** ppExpCntElems
Definition: schurdataman.h:119
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
void BeforePredict(VectorHandler &X, VectorHandler &XP, VectorHandler &XPrev, VectorHandler &XPPrev) const
int mbdyn_METIS_PartGraph(int iTotVertices, int *pXadj, int *pAdjncy, int *pVertexWgts, int *pCommWgts, int *pEdgeWgts, int DataCommSize, int *pParAmgProcs)
Definition: metiswrap.cc:40
int * pAdjncy
Definition: schurdataman.cc:60
#define MYSLEEP(t)
Definition: mysleep.h:49
virtual void AssRes(VectorHandler &ResHdl, doublereal dCoef)
Definition: elman.cc:498
Node ** SearchNode(Node **ppFirst, int dim, unsigned int &label)
#define NO_OP
Definition: myassert.h:74
virtual Elem::Type GetElemType(void) const =0
Elem ** ppMyIntElems
Definition: schurdataman.h:71
VecIter< Elem * > MyElemIter
Definition: schurdataman.h:68
Definition: dataman4.cc:96
void AssRes(VectorHandler &ResHdl, doublereal dCoef)
Definition: schurdataman.cc:93
void AfterPredict(void) const
Node ** ppExpCntNodes
Definition: schurdataman.h:118
integer iTotDofs
Definition: dataman.h:809
std::ostream & Partition(void) const
Definition: output.h:520
void AssJac(MatrixHandler &JacHdl, doublereal dCoef)
Definition: schurdataman.cc:99
void Output(bool force=false) const
doublereal dRestartTime
Definition: dataman.h:166
virtual const InducedVelocity * pGetInducedVelocity(void) const
Definition: aerodyn.cc:791
DriveCaller * pOutputMeter
Definition: dataman.h:179
unsigned int * pLabelsList
Definition: schurdataman.h:84
virtual void MakeRestart(void)
Definition: dataman.cc:699
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
doublereal dLastRestartTime
Definition: dataman.h:173
int * pXadj
Definition: schurdataman.cc:59
int GetDescription(void)
Definition: parser.cc:730
integer * pMyIntDofs
Definition: schurdataman.h:83
struct DataManager::NodeDataStructure NodeData[Node::LASTNODETYPE]
Type
Definition: node.h:71
#define DEBUGCOUT(msg)
Definition: myassert.h:232
virtual void AssJac(MatrixHandler &JacHdl, doublereal dCoef)
Definition: elman.cc:392
VectorHandler * pXPrimeCurr
Definition: dataman.h:109
integer iRestartIterations
Definition: dataman.h:165
Node ** ppIntNodes
Definition: schurdataman.h:80
#define ASSERT(expression)
Definition: colamd.c:977
Elem ** ppFindElem(Elem::Type Typ, unsigned int uElem) const
Definition: elman.cc:629
KeyWords
Definition: dataman4.cc:94
virtual void BeforePredict(VectorHandler &, VectorHandler &, VectorHandler &, VectorHandler &) const
Definition: simentity.cc:82
virtual unsigned int iGetNumDof(void) const
Definition: elem.cc:118
float MAX
Definition: ann_in.c:96
integer iCurrRestartIter
Definition: dataman.h:172
void OutputPartition(void)
unsigned int iNum
Definition: dataman.h:618
Definition: solver.h:78
virtual doublereal dGet(const doublereal &dVar) const =0
virtual void GetConnectedNodes(std::vector< const Node * > &connectedNodes) const
Definition: elem.h:243
virtual unsigned int iGetNumDof(void) const =0
void Init(const T *p, unsigned i)
Definition: veciter.h:79
const char * psNodeNames[]
Definition: enums.cc:372
Dof * pGetDofsList(void) const
virtual bool IsArg(void)
Definition: parser.cc:807
Definition: elem.h:75
virtual Node::Type GetNodeType(void) const =0
Type
Definition: elem.h:91
#define SAFENEWARRNOFILL(pnt, item, sz)
Definition: mynewmem.h:704
virtual void Output(OutputHandler &OH) const
Definition: output.cc:870
integer HowManyDofs(DofType who) const
const char * psElemNames[]
Definition: enums.cc:39
SchurDataManager(MBDynParser &HP, unsigned OF, Solver *pS, doublereal dInitialTime, const char *sOutputFileName, const char *sInputFileName, bool bAbortAfterInput)
Definition: schurdataman.cc:70
virtual int GetWord(void)
Definition: parser.cc:1083
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: simentity.cc:120
void Pack(int *pList, int dim)
void chaco_interface(const int iTotVertices, int *start, int *adjacency, int *vertex_weights, int *comm_weights, int *edge_weights, const int num_processors, int *pParAmgProcs)
Definition: chacowrap.cc:63
integer * pLocalDofs
Definition: schurdataman.h:76
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
virtual void AfterPredict(VectorHandler &X, VectorHandler &XP)
Definition: simentity.cc:91
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
void AfterConvergence(void) const
bool PartitionOpen(void)
Definition: output.cc:630
Definition: dofown.h:59
unsigned int GetLabel(void) const
Definition: withlab.cc:62
struct DataManager::ElemDataStructure ElemData[Elem::LASTELEMTYPE]
Node * pFindNode(Node::Type Typ, unsigned int uNode) const
Definition: nodeman.cc:179
void CreatePartition(void)
const integer iDefaultMaxNodesPerElem
Definition: schurdataman.cc:65
const integer iDefaultMaxConnectionsPerVertex
Definition: schurdataman.cc:64
VariableSubMatrixHandler * pWorkMat
Definition: dataman.h:632
void InitList(int *list, int dim, int value)
VectorHandler * pXCurr
Definition: dataman.h:108
virtual void Update(const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: simentity.cc:98
const int ADJ_UNDEFINED
Definition: schurdataman.cc:67