35 ss <<
"Horror, horror! Trying to create a matrix with size " <<
nCompt_
37 dump(ss.str(),
"ERROR");
38 throw range_error(
"Expected greater than 0.");
45 for (
unsigned int i = 0; i <
nCompt_; i++ )
46 Ga_.push_back( 2.0 / tree[ i ].Ra );
74 const vector< unsigned int >& A,
75 const vector< unsigned int >& B )
77 if ( A.empty() || B.empty() )
80 return A[ 0 ] < B[ 0 ];
87 for (
unsigned int i = 0; i <
nCompt_; ++i )
89 const vector< unsigned int >& c = ( *tree_ )[ i ].children;
96 int diff = ( int )( c[ 0 ] ) - i;
98 if ( diff == 1 || diff == -1 )
109 vector< vector< unsigned int > >::iterator
group;
111 sort( group->begin(), group->end() );
120 for (
unsigned int c = 0; c < group->size() - 1; ++c )
122 index = ( *group )[ c ];
123 rank = group->size() - c - 1;
135 const vector< TreeNodeStruct >& node = *
tree_;
139 for (
unsigned int i = 0; i <
nCompt_; ++i )
141 node[ i ].Cm / (
dt_ / 2.0 ) +
145 vector< JunctionStruct >::iterator junction =
junction_.begin();
146 for (
unsigned int i = 0; i < nCompt_ - 1; ++i )
150 i == junction->index )
158 gij = gi * gj / ( gi + gj );
160 HS_[ 4 * i + 1 ] = -gij;
161 HS_[ 4 * i + 2 ] += gij;
162 HS_[ 4 * i + 6 ] += gij;
165 vector< vector< unsigned int > >::iterator
group;
166 vector< unsigned int >::iterator i;
171 for ( i = group->begin(); i != group->end(); ++i )
174 for ( i = group->begin(); i != group->end(); ++i )
178 HS_[ 4 * *i + 2 ] += gi * ( 1.0 - gi / gsum );
183 vector< unsigned int >::iterator j;
184 unsigned int size = 0;
188 rank = group->size() - 1;
189 size += rank * ( rank + 1 );
198 for ( i = group->begin(); i != group->end(); ++i )
201 for ( i = group->begin(); i != group->end() - 1; ++i )
203 int base =
HJ_.size();
205 for ( j = i + 1; j != group->end(); ++j )
207 gij =
Ga_[ *i ] *
Ga_[ *j ] / gsum;
209 HJ_.push_back( -gij );
210 HJ_.push_back( -gij );
219 for (
unsigned int i = 0; i <
nCompt_; ++i )
220 HS_[ 4 * i ] =
HS_[ 4 * i + 2 ];
230 unsigned int farIndex;
232 vector< JunctionStruct >::iterator junction;
241 index = junction->index;
242 rank = junction->rank;
246 const vector< unsigned int >&
group =
254 farIndex = group[ group.size() - 1 ];
258 else if ( rank == 2 )
263 farIndex = group[ group.size() - 2 ];
268 farIndex = group[ group.size() - 1 ];
275 int start = group.size() - rank;
276 for (
unsigned int j = 0; j < rank; ++j )
278 farIndex = group[ start + j ];
283 operand_.push_back( base + 2 * j + 1 );
286 operand_.push_back(
HS_.begin() + 4 * farIndex + 3 );
288 operand_.push_back( base + 2 * j + 1 );
298 target = base + 2 * rank;
299 for (
unsigned int i = 1; i < rank; ++i )
301 above = base + 2 * i;
302 for (
unsigned int j = 0; j < rank - i; ++j )
315 target = base + 2 * rank + 1;
317 for (
unsigned int i = 1; i < rank; ++i )
319 left = base + 2 * i + 1;
320 for (
unsigned int j = 0; j < rank - i; ++j )
331 if ( i == rank - 1 && j == rank - i - 1 )
345 if ( junction->rank < 3 )
348 index = junction->index;
349 rank = junction->rank;
353 const vector< unsigned int >&
group =
356 unsigned int start = group.size() - rank;
357 for (
unsigned int j = 0; j < rank; ++j )
359 farIndex = group[ start + j ];
388 return HS_[ 4 * row ];
390 unsigned int smaller = row < col ? row : col;
391 unsigned int bigger = row > col ? row : col;
395 if ( bigger - smaller == 1 )
396 return HS_[ 4 * smaller + 1 ];
404 unsigned int groupNumber =
groupNumber_.find( smaller )->second;
405 const vector< unsigned int >&
group =
coupled_[ groupNumber ];
406 unsigned int location, size;
407 unsigned int smallRank, bigRank;
409 if ( find( group.begin(), group.end(), bigger ) != group.end() )
412 for (
int i = 0; i < static_cast< int >( groupNumber ); ++i )
415 location += size * ( size - 1 );
419 smallRank = group.end() - find( group.begin(), group.end(), smaller ) - 1;
420 bigRank = group.end() - find( group.begin(), group.end(), bigger ) - 1;
421 location += size * ( size - 1 ) - smallRank * ( smallRank + 1 );
422 location += 2 * ( smallRank - bigRank - 1 );
424 if ( row == smaller )
425 return HJ_[ location ];
427 return HJ_[ location + 1 ];
438 return HS_[ 4 * row + 3 ];
451 unsigned int size = m.
getSize();
454 for (
unsigned int i = 0; i < size; i++ )
456 for (
unsigned int j = 0; j < size; j++ )
457 s << setw( 12 ) << setprecision( 5 ) << m.
getA( i, j );
462 for (
unsigned int i = 0; i < size; i++ )
466 for (
unsigned int i = 0; i < size; i++ )
467 s << m.
getB( i ) <<
"\n";
478 void testHinesMatrix()
480 vector< int* > childArray;
481 vector< unsigned int > childArraySize;
503 int childArray_1[ ] =
527 childArray.push_back( childArray_1 );
528 childArraySize.push_back(
sizeof( childArray_1 ) /
sizeof(
int ) );
542 int childArray_2[ ] =
550 childArray.push_back( childArray_2 );
551 childArraySize.push_back(
sizeof( childArray_2 ) /
sizeof(
int ) );
565 int childArray_3[ ] =
573 childArray.push_back( childArray_3 );
574 childArraySize.push_back(
sizeof( childArray_3 ) /
sizeof(
int ) );
588 int childArray_4[ ] =
596 childArray.push_back( childArray_4 );
597 childArraySize.push_back(
sizeof( childArray_4 ) /
sizeof(
int ) );
612 int childArray_5[ ] =
622 childArray.push_back( childArray_5 );
623 childArraySize.push_back(
sizeof( childArray_5 ) /
sizeof(
int ) );
638 int childArray_6[ ] =
649 childArray.push_back( childArray_6 );
650 childArraySize.push_back(
sizeof( childArray_6 ) /
sizeof(
int ) );
656 int childArray_7[ ] =
661 childArray.push_back( childArray_7 );
662 childArraySize.push_back(
sizeof( childArray_7 ) /
sizeof(
int ) );
668 int childArray_8[ ] =
675 childArray.push_back( childArray_8 );
676 childArraySize.push_back(
sizeof( childArray_8 ) /
sizeof(
int ) );
682 int childArray_9[ ] =
689 childArray.push_back( childArray_9 );
690 childArraySize.push_back(
sizeof( childArray_9 ) /
sizeof(
int ) );
696 vector< TreeNodeStruct > tree;
703 vector< vector< double > > matrix;
705 double epsilon = 1e-17;
710 unsigned int arraySize;
711 for (
unsigned int cell = 0; cell < childArray.size(); cell++ )
713 array = childArray[ cell ];
714 arraySize = childArraySize[ cell ];
715 nCompt = count( array, array + arraySize, -1 );
719 tree.resize( nCompt );
720 for ( i = 0; i < nCompt; ++i )
722 tree[ i ].Ra = 15.0 + 3.0 * i;
723 tree[ i ].Rm = 45.0 + 15.0 * i;
724 tree[ i ].Cm = 500.0 + 200.0 * i * i;
728 for (
unsigned int a = 0; a < arraySize; a++ )
729 if ( array[ a ] == -1 )
732 tree[ count ].children.push_back( array[ a ] );
741 for ( i = 0; i < nCompt; ++i )
742 for ( j = 0; j < nCompt; ++j )
745 error <<
"Testing Hines' Matrix: Cell# "
746 << cell + 1 <<
", entry (" << i <<
", " << j <<
")";
748 fabs( matrix[ i ][ j ] - H.
getA( i, j ) ) < epsilon,
753 cout <<
"." << flush;
757 #endif // DO_UNIT_TESTS
double getVMid(unsigned int row) const
#define ASSERT(unused, message)
void makeFullMatrix(const vector< TreeNodeStruct > &tree, double dt, vector< vector< double > > &matrix)
static SrcFinfo0 * group()
map< unsigned int, unsigned int > groupNumber_
vector< double > VMid_
middle of a time step.
void setup(const vector< TreeNodeStruct > &tree, double dt)
bool groupCompare(const vector< unsigned int > &A, const vector< unsigned int > &B)
double getB(unsigned int row) const
vector< vdIterator > backOperand_
int stage_
reached. Used in getA.
vector< vdIterator > operand_
const vector< TreeNodeStruct > * tree_
setup.
vector< vector< unsigned int > > coupled_
map< unsigned int, vdIterator > operandBase_
ostream & operator<<(ostream &s, const HinesMatrix &m)
void makeOperands()
elimination easier.
unsigned int getSize() const
vector< JunctionStruct > junction_
vector< double >::iterator vdIterator
double getA(unsigned int row, unsigned int col) const