ZeroDimensionalCircuitData.cpp 40.8 KB
Newer Older
1
2
//@HEADER
/*
Cristiano Malossi's avatar
Cristiano Malossi committed
3
*******************************************************************************
4

Cristiano Malossi's avatar
Cristiano Malossi committed
5
6
    Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
    Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7

Cristiano Malossi's avatar
Cristiano Malossi committed
8
    This file is part of LifeV.
9

Cristiano Malossi's avatar
Cristiano Malossi committed
10
11
12
13
    LifeV is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
14

Cristiano Malossi's avatar
Cristiano Malossi committed
15
16
17
18
    LifeV is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.
19

Cristiano Malossi's avatar
Cristiano Malossi committed
20
21
    You should have received a copy of the GNU Lesser General Public License
    along with LifeV.  If not, see <http://www.gnu.org/licenses/>.
22

Cristiano Malossi's avatar
Cristiano Malossi committed
23
24
*******************************************************************************
*/
25
26
27
28
//@HEADER

/*!
 *  @file
Cristiano Malossi's avatar
Cristiano Malossi committed
29
 *  @brief File containing a class for 0D model circuit data handling.
30
 *
Cristiano Malossi's avatar
Cristiano Malossi committed
31
32
 *  @date 26-09-2011
 *  @author Mahmoud Jafargholi <mahmoud.jafargholi@epfl.ch>
33
 *
Cristiano Malossi's avatar
Cristiano Malossi committed
34
35
 *  @contributors Cristiano Malossi <cristiano.malossi@epfl.ch>
 *  @mantainer    Cristiano Malossi <cristiano.malossi@epfl.ch>
36
37
 */

38
#include <lifev/zero_dimensional/solver/ZeroDimensionalCircuitData.hpp>
39

40
namespace LifeV {
41
42

// ===================================================
43
// Methods
44
// ===================================================
45
void ZeroDimensionalElement::showMe( const Int& /*flag*/)
46
{
47
48
49
    cout << "Id = " << id();
    cout << "\t type = ";
    cout << enum2string( type() );
50
51
}

52
const std::string ZeroDimensionalElement::enum2string( const ZeroDimensionalElementType & type )
53
54
{
    std::string output;
55
56
57
    switch ( type )
    {
        case resistor:
58
        {
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
            output = "resistor";
            break;
        }
        case capacitor:
        {
            output = "capacitor";
            break;
        }
        case inductor:
        {
            output = "inductor";
            break;
        }
        case diode:
        {
            output = "diode";
            break;
        }
        case voltageSource:
        {
            output = "voltageSource";
            break;
        }
        case currentSource:
        {
            output = "currentSource";
            break;
86
        }
87
88
    }
    return ( output );
89
90
91
92
93
}

// ===================================================
// Constructors
// ===================================================
94
95
ZeroDimensionalElementPassive::ZeroDimensionalElementPassive() :
    M_parameter(), M_nodeIndex()
96
{
97
    M_nodeIndex.reserve( 2 );
98
99
}

100
void ZeroDimensionalElementPassive::showMe( const Int& flag )
101
{
102
103
    ZeroDimensionalElement::showMe( flag );
    cout << "\t Node1= " << nodeIndex( 0 ) << "\t Node2= " << nodeIndex( 1 );
104
105
}

106
void ZeroDimensionalElementPassive::connectElement( zeroDimensionalNodeSPtr_Type & Nodes )
107
{
108
109
    Int Index = M_nodeIndex.at( 0 );
    zeroDimensionalNodePtr_Type theNode = Nodes->nodeListAt( Index );
110

111
    theNode->setElementListIndex( M_id );
112
    Int otherEndIndex = M_nodeIndex.at( 1 );
113
    theNode->setNodeListIndex( otherEndIndex );
114

115
116
    Index = M_nodeIndex.at( 1 );
    theNode = Nodes->nodeListAt( Index );
117

118
    theNode->setElementListIndex( M_id );
119
    otherEndIndex = M_nodeIndex.at( 0 );
120
    theNode->setNodeListIndex( otherEndIndex );
121
122
}

123
124
125
126
Real ZeroDimensionalElementPassive::direction( const Int & nodeId ) const
{
    if ( M_nodeIndex.at( 0 ) == nodeId )
    {
127
        return -1.0;
128
129
130
131
    }
    else
    {
        return 1.0;
132
133
134
135
136
137
138
    }
}
// ===================================================
// Constructors
// ===================================================
ZeroDimensionalElementPassiveResistor::ZeroDimensionalElementPassiveResistor()
{
139
    M_type = resistor;
140
141
}

142
void ZeroDimensionalElementPassiveResistor::showMe( const Int& flag )
143
{
144
145
    ZeroDimensionalElementPassive::showMe( flag );
    cout << "\t parameter = " << parameter() << std::endl;
146
147
}

148
149
150
151
void ZeroDimensionalElementPassiveResistor::buildABC( matrix_Type& /*A*/,
                                                      matrix_Type& B,
                                                      vector_Type& C,
                                                      const zeroDimensionalNodeSPtr_Type& Nodes )
152
{
153
154
155
156
157
158
159
160
161
    for ( Int i = 0; i < 2; i++ )
    {
        const Int& theNodeCounter = ( i ) % 2;
        const Int& theOtherNodeCounter = ( i + 1 ) % 2;
        const ZeroDimensionalNode& theNodeTest = *( Nodes->nodeListAt( M_nodeIndex[theNodeCounter] ) );
        const ZeroDimensionalNode& theOtherNodeTest = *( Nodes->nodeListAt( M_nodeIndex[theOtherNodeCounter] ) );
        if ( theNodeTest.type() == unknownNode )
        {
            const ZeroDimensionalNodeUnknown& theNode = *( Nodes->unknownNodeMapAt( M_nodeIndex[theNodeCounter] ) );
162
            const Int& equationRow = theNode.equationRow();
163
164
165
166
167
168
169
170
171
172
173
174
175
176
            B.addToCoefficient( equationRow,
                                theNode.variableIndex(),
                                -M_parameter );
            if ( theOtherNodeTest.type() == unknownNode )
            {
                const ZeroDimensionalNodeUnknown& theOtherNode = *( Nodes->unknownNodeMapAt( M_nodeIndex[theOtherNodeCounter] ) );
                B.addToCoefficient( equationRow,
                                    theOtherNode.variableIndex(),
                                    M_parameter );
            }
            else
            {
                const ZeroDimensionalNodeKnown & theOtherNode = *( Nodes->knownNodeMapAt( M_nodeIndex[theOtherNodeCounter] ) );
                C[equationRow] += ( M_parameter * theOtherNode.voltage() );
177
178
179
180
181
182

            }
        }
    }
}

183
void ZeroDimensionalElementPassiveResistor::extractSolution( const ZeroDimensionalNodeS& Nodes )
184
{
185
186
187
    M_current = ( Nodes.nodeListAt( M_nodeIndex.at( 0 ) )->voltage() - Nodes.nodeListAt( M_nodeIndex.at( 1 ) )->voltage() ) * M_parameter;
    M_deltaCurrent = ( Nodes.nodeListAt( M_nodeIndex.at( 0 ) )->deltaVoltage() - Nodes.nodeListAt( M_nodeIndex.at( 1 ) )->deltaVoltage() )
                    * M_parameter;
188
189
190
191
}
// ===================================================
// Constructors
// ===================================================
192
193
ZeroDimensionalElementPassiveDiode::ZeroDimensionalElementPassiveDiode() :
    M_alpha(), M_beta(), M_forwardBias()
194
{
195
    M_type = diode;
196
197
}

198
void ZeroDimensionalElementPassiveDiode::showMe( const Int& flag )
199
{
200
201
    ZeroDimensionalElementPassiveResistor::showMe( flag );
    cout << "\t FB = " << forwardBias() << "\t alpha = " << alpha() << "\t beta = " << beta() << std::endl;
202
}
203
204
void ZeroDimensionalElementPassiveDiode::calculateEffectiveResistance( const Real& voltage )
{
205

Cristiano Malossi's avatar
Cristiano Malossi committed
206
207
    Real e0 = ( M_beta * exp( M_alpha * ( - M_forwardBias ) ) );
    Real e1 = M_alpha * e0;
208
209
    if ( voltage == 0 )
    {
Cristiano Malossi's avatar
Cristiano Malossi committed
210
        M_parameter = e1;
211
212
213
214
215
    }
    else
    {
        Real current = ( M_beta * exp( M_alpha * ( voltage - M_forwardBias ) ) ) - e0;
        M_parameter = abs( current / voltage );
216
217
    }
}
218
void ZeroDimensionalElementPassiveDiode::extractSolution( const ZeroDimensionalNodeS& Nodes )
219
{
220
221
222
223
224
    Real deltaVoltage = ( Nodes.nodeListAt( M_nodeIndex.at( 0 ) )->voltage() - Nodes.nodeListAt( M_nodeIndex.at( 1 ) )->voltage() );
    calculateEffectiveResistance( deltaVoltage );
    M_current = deltaVoltage * M_parameter;
    M_deltaCurrent = ( Nodes.nodeListAt( M_nodeIndex.at( 0 ) )->deltaVoltage() - Nodes.nodeListAt( M_nodeIndex.at( 1 ) )->deltaVoltage() )
                    * M_parameter;
225
226
}

227
228
229
230
void ZeroDimensionalElementPassiveDiode::buildABC( matrix_Type& A,
                                                   matrix_Type& B,
                                                   vector_Type& C,
                                                   const zeroDimensionalNodeSPtr_Type& Nodes )
231
232
{

233
234
235
236
    const Int& theNodeCounter = 0;
    const Int& theOtherNodeCounter = 1;
    const ZeroDimensionalNode& theNodeTest = *( Nodes->nodeListAt( M_nodeIndex[theNodeCounter] ) );
    const ZeroDimensionalNode& theOtherNodeTest = *( Nodes->nodeListAt( M_nodeIndex[theOtherNodeCounter] ) );
237
    Real deltaVoltage = theNodeTest.voltage() - theOtherNodeTest.voltage();
238
239
    calculateEffectiveResistance( deltaVoltage );
    ZeroDimensionalElementPassiveResistor::buildABC( A, B, C, Nodes );
240
241
242
243
244
245
246
}

// ===================================================
// Constructors
// ===================================================
ZeroDimensionalElementPassiveCapacitor::ZeroDimensionalElementPassiveCapacitor()
{
247
    M_type = capacitor;
248
249
}

250
void ZeroDimensionalElementPassiveCapacitor::showMe( const Int& flag )
251
{
252
253
    ZeroDimensionalElementPassive::showMe( flag );
    cout << "\t parameter = " << parameter() << std::endl;
254
255
}

256
void ZeroDimensionalElementPassiveCapacitor::extractSolution( const ZeroDimensionalNodeS& Nodes )
257
{
258
    M_current = ( Nodes.nodeListAt( M_nodeIndex.at( 0 ) )->deltaVoltage() - Nodes.nodeListAt( M_nodeIndex.at( 1 ) )->deltaVoltage() ) * M_parameter;
259
260
}

261
262
263
264
void ZeroDimensionalElementPassiveCapacitor::buildABC( matrix_Type& A,
                                                       matrix_Type& /*B*/,
                                                       vector_Type& C,
                                                       const zeroDimensionalNodeSPtr_Type& Nodes )
265
{
266
267
268
269
270
271
272
273
274
    for ( Int i = 0; i < 2; i++ )
    {
        const Int& theNodeCounter = ( i ) % 2;
        const Int& theOtherNodeCounter = ( i + 1 ) % 2;
        const ZeroDimensionalNode& theNodeTest = *( Nodes->nodeListAt( M_nodeIndex[theNodeCounter] ) );
        const ZeroDimensionalNode& theOtherNodeTest = *( Nodes->nodeListAt( M_nodeIndex[theOtherNodeCounter] ) );
        if ( theNodeTest.type() == unknownNode )
        {
            const ZeroDimensionalNodeUnknown& theNode = *( Nodes->unknownNodeMapAt( M_nodeIndex[theNodeCounter] ) );
275
            const Int& equationRow = theNode.equationRow();
276
277
278
279
280
281
282
283
284
285
286
287
288
289
            A.addToCoefficient( equationRow,
                                theNode.variableIndex(),
                                -M_parameter );
            if ( theOtherNodeTest.type() == unknownNode )
            {
                const ZeroDimensionalNodeUnknown& theOtherNode = *( Nodes->unknownNodeMapAt( M_nodeIndex[theOtherNodeCounter] ) );
                A.addToCoefficient( equationRow,
                                    theOtherNode.variableIndex(),
                                    M_parameter );
            }
            else
            {
                const ZeroDimensionalNodeKnown & theOtherNode = *( Nodes->knownNodeMapAt( M_nodeIndex[theOtherNodeCounter] ) );
                C[equationRow] += ( M_parameter * theOtherNode.deltaVoltage() );
290
291
292
293
294
295
296
            }
        }
    }
}
// ===================================================
// Constructors
// ===================================================
297
298
ZeroDimensionalElementPassiveInductor::ZeroDimensionalElementPassiveInductor() :
    M_equationRow(), M_variableIndex()
299
{
300
    M_type = inductor;
301
302
}

303
void ZeroDimensionalElementPassiveInductor::showMe( const Int& flag )
304
{
305
306
307
    ZeroDimensionalElementPassive::showMe( flag );
    if ( flag == 0 )
        cout << "\t parameter = " << parameter() << std::endl;
308
    else
309
        cout << "\t EquationRow= " << M_equationRow << "\t VariableIndex= " << M_variableIndex;
310
311
}

312
void ZeroDimensionalElementPassiveInductor::assignVariableIndex( const Int & index )
313
{
314
315
    M_equationRow = index;
    M_variableIndex = index;
316
317

}
318
319
320
321
void ZeroDimensionalElementPassiveInductor::buildABC( matrix_Type& A,
                                                      matrix_Type& B,
                                                      vector_Type& C,
                                                      const zeroDimensionalNodeSPtr_Type& Nodes )
322
{
323
324
325
326
327
328
329
    for ( Int i = 0; i < 2; i++ )
    {
        const Int& theNodeCounter = ( i ) % 2;
        const ZeroDimensionalNode& theNodeTest = *( Nodes->nodeListAt( M_nodeIndex[theNodeCounter] ) );
        if ( theNodeTest.type() == unknownNode )
        {
            const ZeroDimensionalNodeUnknown& theNode = *( Nodes->unknownNodeMapAt( M_nodeIndex[theNodeCounter] ) );
330
            const Int& equationRow = theNode.equationRow();
331
332
333
            B.addToCoefficient( equationRow,
                                M_variableIndex,
                                direction( theNode.id() ) );
334
335
        }
    }
336
337
338
339
    // Write the equations for an inductor
    A.addToCoefficient( M_equationRow,
                        M_variableIndex,
                        1.0 );
340
    {
341
        Int i = 0;
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
        const Int& theNodeCounter = ( i ) % 2;
        const Int& theOtherNodeCounter = ( i + 1 ) % 2;
        const ZeroDimensionalNode& theNodeTest = *( Nodes->nodeListAt( M_nodeIndex[theNodeCounter] ) );
        const ZeroDimensionalNode& theOtherNodeTest = *( Nodes->nodeListAt( M_nodeIndex[theOtherNodeCounter] ) );
        if ( theNodeTest.type() == unknownNode )
        {
            const ZeroDimensionalNodeUnknown& theNode = *( Nodes->unknownNodeMapAt( M_nodeIndex[theNodeCounter] ) );
            B.addToCoefficient( M_equationRow,
                                theNode.variableIndex(),
                                -M_parameter );
        }
        else
        {
            const ZeroDimensionalNodeKnown & theNode = *( Nodes->knownNodeMapAt( M_nodeIndex[theNodeCounter] ) );
            C[M_equationRow] += ( -M_parameter * theNode.voltage() );
357
358
        }

359
360
361
362
363
364
365
366
367
368
369
        if ( theOtherNodeTest.type() == unknownNode )
        {
            const ZeroDimensionalNodeUnknown& theOtherNode = *( Nodes->unknownNodeMapAt( M_nodeIndex[theOtherNodeCounter] ) );
            B.addToCoefficient( M_equationRow,
                                theOtherNode.variableIndex(),
                                M_parameter );
        }
        else
        {
            const ZeroDimensionalNodeKnown& theOtherNode = *( Nodes->knownNodeMapAt( M_nodeIndex[theOtherNodeCounter] ) );
            C[M_equationRow] += ( +M_parameter * theOtherNode.voltage() );
370
371
372
373
374
375
376
377
        }
    }

}

// ==================================================
// Constructors
// ===================================================
378
379
ZeroDimensionalElementSource::ZeroDimensionalElementSource() :
    M_nodeIndex()
380
381
382
383

{
}

384
void ZeroDimensionalElementSource::showMe( const Int& flag )
385
{
386
387
    ZeroDimensionalElement::showMe( flag );
    cout << "\t Node1= " << nodeIndex() << std::endl;
388
389
390
391
392
}

// ===================================================
// Constructors
// ===================================================
393
394
ZeroDimensionalElementVoltageSource::ZeroDimensionalElementVoltageSource() :
    M_voltage(), M_deltaVoltage()
395
{
396
    M_type = voltageSource;
397
398
399

}

400
void ZeroDimensionalElementVoltageSource::connectElement( zeroDimensionalNodeSPtr_Type & Nodes )
401
{
402
    zeroDimensionalNodeKnownPtr_Type theNode = Nodes->knownNodeMapAt( M_nodeIndex );
403
    theNode->setElementListIndex( M_id );
404
    Int otherEndIndex = -1;// In fact this element is not connecting this node to a new node, only filling the vector
405
    theNode->setNodeListIndex( otherEndIndex );
406
}
407
408
409
410
void ZeroDimensionalElementVoltageSource::calculateCurrent( const ZeroDimensionalNodeS& Nodes,
                                                            const ZeroDimensionalElementS& Elements )
{
    ZeroDimensionalNodeKnown& theNode = *Nodes.knownNodeMapAt( M_nodeIndex );
411
    const vecInt_Type& indexList = theNode.elementListIndex();
412
    Real current = 0.0;
413
    Int length = indexList.size();
414
415
416
417
418
419
420
    for ( Int i = 0; i < length; ++i )
    {
        const ZeroDimensionalElement& theElement = *Elements.elementListAt( indexList.at( i ) );
        if ( theElement.id() != M_id )
        {
            Real tmp = theElement.current() * theElement.direction( theNode.id() );
            current += tmp;
421
422
423
424
425
426
427
428
429
        }
    }
    M_current = current;
}
// ===================================================
// Constructors
// ===================================================
ZeroDimensionalElementCurrentSource::ZeroDimensionalElementCurrentSource()
{
430
    M_type = currentSource;
431
432
}

433
void ZeroDimensionalElementCurrentSource::connectElement( zeroDimensionalNodeSPtr_Type & Nodes )
434
{
435
436
    Int Index = M_nodeIndex;
    zeroDimensionalNodePtr_Type theNode = Nodes->nodeListAt( Index );
437
    theNode->setElementListIndex( M_id );
438
    Int otherEndIndex = -1;// In fact this element is not connecting this node to a new node, only filling the vector
439
    theNode->setNodeListIndex( otherEndIndex );
440
441
}

442
443
444
445
void ZeroDimensionalElementCurrentSource::buildABC( matrix_Type& /*A*/,
                                                    matrix_Type& /*B*/,
                                                    vector_Type& C,
                                                    const zeroDimensionalNodeSPtr_Type& Nodes )
446
{
447
448
449
450
    const ZeroDimensionalNode& theNodeTest = *( Nodes->nodeListAt( M_nodeIndex ) );
    if ( theNodeTest.type() == unknownNode )
    {
        const ZeroDimensionalNodeUnknown& theNode = *( Nodes->unknownNodeMapAt( M_nodeIndex ) );
451
        const Int& equationRow = theNode.equationRow();
452
453
454
455
456
457
        C[equationRow] += ( -M_current );
    }
    else
    {
        std::cerr << "Error at ZeroDimensionalElementCurrentSource::buildABC, source connected to voltage";
        std::exit( -1 );
458
459
460
461
462
463
    }
}

// ===================================================
// Constructors
// ===================================================
464
465
ZeroDimensionalNode::ZeroDimensionalNode() :
    M_id(), M_type(), M_currentBalance( 0 ), M_elementListIndex(), M_nodeListIndex(), M_voltage(), M_deltaVoltage()
466
467
468
{
}

469
void ZeroDimensionalNode::showMe( const Int& flag )
470
{
471
472
473
474
475
476
477
478
479
    cout << "Id = " << id();
    cout << "\t type = ";
    cout << enum2string( type() ) << "\t";
    if ( flag == 0 )
    {
        for ( UInt i = 0; i < M_elementListIndex.size(); i++ )
            cout << "{" << M_elementListIndex.at( i ) << "," << M_nodeListIndex.at( i ) << "}  ";
        cout << std::endl;
    }
480
481
}

482
const std::string ZeroDimensionalNode::enum2string( const ZeroDimensionalNodeType & type ) const
483
484
{
    std::string output;
485
486
487
    switch ( type )
    {
        case knownNode:
488
        {
489
490
            output = "knownNode";
            break;
491
        }
492
493
494
495
496
497
498
499
        case unknownNode:
        {
            output = "unknownNode";
            break;
        }

    }
    return ( output );
500
501
}

502
503
504
void ZeroDimensionalNode::calculateCurrentBalance( const ZeroDimensionalElementS& Elements )
{
    Real currentBalance = 0.0;
505
    Int length = M_elementListIndex.size();
506
507
508
509
    for ( Int i = 0; i < length; ++i )
    {
        const ZeroDimensionalElement& theElement = *Elements.elementListAt( M_elementListIndex.at( i ) );
        currentBalance += theElement.current() * theElement.direction( M_id );
510
511
512
513
514
515
516
    }
    M_currentBalance = currentBalance;
}

// ===================================================
// Constructors
// ===================================================
517
518
ZeroDimensionalNodeUnknown::ZeroDimensionalNodeUnknown() :
    M_variableIndex(), M_equationRow()
519
520
521
{
    M_type = unknownNode;
}
522
void ZeroDimensionalNodeUnknown::assignVariableIndex( const Int & index )
523
{
524
525
    M_equationRow = index;
    M_variableIndex = index;
526
527
528

}

529
void ZeroDimensionalNodeUnknown::showMe( const Int& flag )
530
{
531
532
533
    ZeroDimensionalNode::showMe( flag );
    if ( flag == 1 )
        cout << "\t EquationRow= " << M_equationRow << "\t VariableIndex= " << M_variableIndex;
534
535
536
537
538
539
}

// ===================================================
// Constructors
// ===================================================
ZeroDimensionalNodeKnown::ZeroDimensionalNodeKnown()
540
541
542
{
    M_type = knownNode;
}
543
544
ZeroDimensionalNodeKnown::ZeroDimensionalNodeKnown( const zeroDimensionalElementVoltageSourcePtr_Type& theElement )
{
545
    M_element = theElement;
546
547
548
549
550
551
    M_type = knownNode;
}

// ===================================================
// Constructors
// ===================================================
552
553
554
555
556
557
558
ZeroDimensionalElementS::ZeroDimensionalElementS() :
    M_elementList( new vecZeroDimensionalElementPtr_Type ), M_resistorList( new vecZeroDimensionalElementPassiveResistorPtr_Type ),
                    M_capacitorList( new vecZeroDimensionalElementPassiveCapacitorPtr_Type ),
                    M_inductorList( new vecZeroDimensionalElementPassiveInductorPtr_Type ),
                    M_diodeList( new vecZeroDimensionalElementPassiveDiodePtr_Type ),
                    M_currentSourceList( new vecZeroDimensionalElementCurrentSourcePtr_Type ),
                    M_voltageSourceList( new vecZeroDimensionalElementVoltageSourcePtr_Type ), M_voltageSourceMap( new mapVoltageSource_Type )
559
560
561
{
}

562
void ZeroDimensionalElementS::showMe( const Int& flag )
563
{
564
    cout << "=============== Show all ZeroDimensional Elements ===========" << std::endl;
565
566
567
    for ( iterZeroDimensionalElement_Type theElement = M_elementList->begin(); theElement != M_elementList->end(); theElement++ )
    {
        ( *theElement )->showMe( flag );
568
569
570
571
572
573
    }
}

// ===================================================
// Constructors
// ===================================================
574
575
576
577
ZeroDimensionalNodeS::ZeroDimensionalNodeS() :
    M_nodeList( new vecZeroDimensionalNodePtr_Type ), M_unknownNodeList( new vecZeroDimensionalNodeUnknownPtr_Type ),
                    M_knownNodeList( new vecZeroDimensionalNodeKnownPtr_Type ), M_knownNodeMap( new mapNodeKnown_Type ),
                    M_unknownNodeMap( new mapNodeUnknown_Type )
578
579
580
{
}

581
void ZeroDimensionalNodeS::showMe( const Int& flag )
582
{
583
    cout << "=============== Show all ZeroDimensional Nodes    ===========" << std::endl;
584
585
586
    for ( iterZeroDimensionalNode_Type theNode = M_nodeList->begin(); theNode != M_nodeList->end(); theNode++ )
    {
        ( *theNode )->showMe( flag );
587
588
589
590
591
592
    }
}

// ===================================================
// Constructors
// ===================================================
593
594
ZeroDimensionalCircuitData::ZeroDimensionalCircuitData() :
    M_Elements( new ZeroDimensionalElementS ), M_Nodes( new ZeroDimensionalNodeS ), M_bc()
595
596
597
598

{
}

599
void ZeroDimensionalCircuitData::showMe( const Int& flag )
600
{
601
602
    M_Elements->showMe( flag );
    M_Nodes->showMe( flag );
603
604
605
606
607
}

// ===================================================
// Methods
// ===================================================
608
609
void ZeroDimensionalCircuitData::buildCircuit( const char *fileName,
                                               bcPtr_Type bc )
610
{
611
    using namespace std;
612
613
614
615
616
    ifstream infile;
    string stringline1;
    string stringline2;
    string stringline3;
    string stringtmp;
617
618
    stringstream stringStreamLine1( stringstream::in | stringstream::out );
    stringstream stringStreamLine2( stringstream::in | stringstream::out );
619
620
621
622
623
624
625
626
627
628
629

    bool boolElementType[ZERO_DIMENTIONAL_DEFINED_ELEMENTS];
    Int numberOfNodes = 0;
    Int numberOfElements = 0;
    Int numberOfTerminalNodes = 0;
    Int ID = 0;
    Int Node1 = 0;
    Int Node2 = 0;
    Real parameter1 = 0.0;
    Real forwardBias = 0.0;
    Real alpha = 0.0;
630
    Real beta = 0.0;
631
632


633
    std::vector< ZeroDimensionalNodeType > nodesType;
634
635
    vecInt_Type nodesConnectingSource;
    vecInt_Type terminalNodes;
636

637
638
639
    infile.open( fileName,
                 ifstream::in );
    if ( infile.is_open() )
640
    {
641
642
        // ----------Read the first line----------------------------
        getline( infile, stringline1 );
643
644
        stringStreamLine1 << stringline1;
        stringStreamLine1 >> stringtmp;
645
        numberOfElements = std::atoi( stringtmp.c_str() );
646
        stringStreamLine1 >> stringtmp;
647
        numberOfNodes = std::atoi( stringtmp.c_str() );
648
        stringStreamLine1 >> stringtmp;
649
650
        numberOfTerminalNodes = std::atoi( stringtmp.c_str() );
        getline( infile, stringline2 );
651
        Int nodeId = -1;
652
653

        // ----------Read the second line---------------------------
654
        stringStreamLine2 << stringline2;
655
656
        for ( Int i = 0; i < numberOfTerminalNodes; i++ )
        {
657
            stringStreamLine2 >> stringtmp;
658
659
            nodeId = std::atoi( stringtmp.c_str() );
            terminalNodes.push_back( nodeId );
660
        }
661
662
663
664
665
        //---------------------------------------------------------
        for ( Int i = 0; i < numberOfNodes; i++ )
        {
            nodesType.push_back( unknownNode );
            nodesConnectingSource.push_back( -1 );
666
        }
667
668
        // ----------Read Elements-----------------------------------
        for ( Int i = 0; i < numberOfElements; i++ )
669
        {
670
671
            stringstream stringStreamLine3( stringstream::in | stringstream::out );
            getline( infile, stringline3 );
672
673
            stringStreamLine3 << stringline3;
            stringStreamLine3 >> stringtmp;
674
            ID = std::atoi( stringtmp.c_str() );
675
            stringStreamLine3 >> stringtmp;
676
677
678
679
680
681
            boolElementType[0] = stringtmp.compare( "resistor" );
            boolElementType[1] = stringtmp.compare( "capacitor" );
            boolElementType[2] = stringtmp.compare( "inductor" );
            boolElementType[3] = stringtmp.compare( "voltageSource" );
            boolElementType[4] = stringtmp.compare( "currentSource" );
            boolElementType[5] = stringtmp.compare( "diode" );
682
683

            stringStreamLine3 >> stringtmp;
684
            Node1 = std::atoi( stringtmp.c_str() );
685

686
            if ( !boolElementType[0] )//resistor
687
688
            {
                stringStreamLine3 >> stringtmp;
689
                Node2 = std::atoi( stringtmp.c_str() );
690
                stringStreamLine3 >> stringtmp;
691
692
                parameter1 = std::atof( stringtmp.c_str() );
                createElementResistor( ID, Node1, Node2, parameter1 );
693
            }
694
            if ( !boolElementType[1] )//capacitor
695
696
            {
                stringStreamLine3 >> stringtmp;
697
                Node2 = std::atoi( stringtmp.c_str() );
698
                stringStreamLine3 >> stringtmp;
699
700
                parameter1 = std::atof( stringtmp.c_str() );
                createElementCapacitor( ID, Node1, Node2, parameter1 );
701
            }
702
            if ( !boolElementType[2] )//inductor
703
704
            {
                stringStreamLine3 >> stringtmp;
705
                Node2 = std::atoi( stringtmp.c_str() );
706
                stringStreamLine3 >> stringtmp;
707
708
                parameter1 = std::atof( stringtmp.c_str() );
                createElementInductor( ID, Node1, Node2, parameter1 );
709
            }
710
            if ( !boolElementType[5] )//diode
711
712
            {
                stringStreamLine3 >> stringtmp;
713
                Node2 = std::atoi( stringtmp.c_str() );
714
                stringStreamLine3 >> stringtmp;
715
                forwardBias = std::atof( stringtmp.c_str() );
716
                stringStreamLine3 >> stringtmp;
717
                alpha = std::atof( stringtmp.c_str() );
718
                stringStreamLine3 >> stringtmp;
719
720
                beta = std::atof( stringtmp.c_str() );
                createElementDiode( ID, Node1, Node2, forwardBias, alpha, beta );
721
722
723
724
725
726
727
            }
        }
        infile.close();
    }
    else
    {
        cerr << "Error opening circuit file";
728
        exit( -1 );
729
730
    }

731
732
733
734
735
736
737
738
739
740
741
742
743
744
    for ( iterVecInt_Type theNodeId = terminalNodes.begin(); theNodeId != terminalNodes.end(); theNodeId++ )
    {
        // create Source Elements -----------------------------------------------
        switch ( bc->bc( *theNodeId ).bcType() )
        {
            case Voltage://create voltage source
                nodesConnectingSource.at( *theNodeId ) = createElementVoltageSource( *theNodeId );
                nodesType.at( *theNodeId ) = knownNode;
                break;
            case Current:// current source
                createElementCurrentSource( *theNodeId );
                break;
            default:
                break;
745
746
747
        }
    }
    //create Nodes --------------------------------------------------------
748
    for ( Int i = 0; i < numberOfNodes; i++ )
749
    {
750
        if ( nodesType.at( i ) == knownNode )
751
        {
752
753
            //Create known Node and connect the voltage source to the node
            createKnownNode( i, M_Elements->voltageSourceMap( nodesConnectingSource.at( i ) ) );
754
755
756
        }
        else
        {
757
            createUnknownNode( i );
758
759
760
        }
    }
    //connect elements to the nodes
761
762
763
    for ( Int i = 0; i < M_Elements->elementCounter(); i++ )
    {
        M_Elements->elementListAt( i )->connectElement( M_Nodes );
764
765
    }
    //set BC to source elements
766
    fixBC( bc );
767
768
}

769
770
771
772
void ZeroDimensionalCircuitData::createElementResistor( Int ID,
                                                        Int node1,
                                                        Int node2,
                                                        Real parameter )
773
{
774
    zeroDimensionalElementPassiveResistorPtr_Type theElement( new ZeroDimensionalElementPassiveResistor() );
775
    theElement->setId( ID );
776
777
778
779
    if ( M_Elements->elementCounter() != ID )
    {
        cerr << "Error: Element Id error at  " << ID;
        exit( -1 );
780
    }
781
782
    theElement->setNodeIndex( node1 );
    theElement->setNodeIndex( node2 );
783
784
785
786
    if ( parameter <= 0 )
    {
        cerr << "Error: Resistance value <=0, ID =  " << ID;
        exit( -1 );
787
    }
788
    theElement->setParameter( 1.0 / parameter );
789
    M_Elements->setelementList( theElement );
790
    M_Elements->setResistorList( theElement );
791
792
793
794
795
796
797
}
void ZeroDimensionalCircuitData::createElementCapacitor( Int ID,
                                                         Int node1,
                                                         Int node2,
                                                         Real parameter )
{
    zeroDimensionalElementPassiveCapacitorPtr_Type theElement( new ZeroDimensionalElementPassiveCapacitor() );
798
    theElement->setId( ID );
799
800
801
802
    if ( M_Elements->elementCounter() != ID )
    {
        cerr << "Error: Element Id error at  " << ID;
        exit( -1 );
803
    }
804
805
806
    theElement->setNodeIndex( node1 );
    theElement->setNodeIndex( node2 );
    theElement->setParameter( parameter );
807
    M_Elements->setelementList( theElement );
808
    M_Elements->setCapacitorList( theElement );
809
810
811
812
813
}
void ZeroDimensionalCircuitData::createElementInductor( Int ID,
                                                        Int node1,
                                                        Int node2,
                                                        Real parameter )
814
{
815
    zeroDimensionalElementPassiveInductorPtr_Type theElement( new ZeroDimensionalElementPassiveInductor() );
816
    theElement->setId( ID );
817

818
819
820
821
    if ( M_Elements->elementCounter() != ID )
    {
        cerr << "Error: Element Id error at  " << ID;
        exit( -1 );
822
    }
823
824
825
    theElement->setNodeIndex( node1 );
    theElement->setNodeIndex( node2 );
    theElement->setParameter( 1.0 / parameter );
826
    M_Elements->setelementList( theElement );
827
    M_Elements->setInductorList( theElement );
828
829
830
831
832
833
834
835
836
}
void ZeroDimensionalCircuitData::createElementDiode( Int ID,
                                                     Int node1,
                                                     Int node2,
                                                     Real forwardBias,
                                                     Real alpha,
                                                     Real beta )
{
    zeroDimensionalElementPassiveDiodePtr_Type theElement( new ZeroDimensionalElementPassiveDiode() );
837
    theElement->setId( ID );
838
839
840
841
    if ( M_Elements->elementCounter() != ID )
    {
        cerr << "Error: Element Id error at  " << ID;
        exit( -1 );
842
    }
843
844
845
    theElement->setNodeIndex( node1 );
    theElement->setNodeIndex( node2 );
    theElement->setParameter( 0 );
846
847
848
849
    theElement->setforwardBias( forwardBias );
    theElement->setalpha( alpha );
    theElement->setbeta( beta );
    M_Elements->setelementList( theElement );
850
    M_Elements->setDiodeList( theElement );
851
852
853
854
855
}

Int ZeroDimensionalCircuitData::createElementVoltageSource( Int node1 )
{
    zeroDimensionalElementVoltageSourcePtr_Type theElement( new ZeroDimensionalElementVoltageSource() );
856
857
    theElement->setId( M_Elements->elementCounter() );
    theElement->setNodeIndex( node1 );
858
    M_Elements->setelementList( theElement );
859
860
    M_Elements->setVoltageSourceList( theElement );
    M_Elements->setVoltageSourceMap( theElement->id(), theElement );
861
862
    return theElement->id();
}
863
void ZeroDimensionalCircuitData::createElementCurrentSource( Int node1 )
864
{
865
    zeroDimensionalElementCurrentSourcePtr_Type theElement( new ZeroDimensionalElementCurrentSource() );
866
867
    theElement->setId( M_Elements->elementCounter() );
    theElement->setNodeIndex( node1 );
868
    M_Elements->setelementList( theElement );
869
    M_Elements->setCurrentSourceList( theElement );
870
871
}

872
void ZeroDimensionalCircuitData::createUnknownNode( const Int & id )
873
{
874
    zeroDimensionalNodeUnknownPtr_Type theNode( new ZeroDimensionalNodeUnknown() );
875
    theNode->setId( id );
876
877
878
    M_Nodes->setnodeList( theNode );
    M_Nodes->setunknownNodeList( theNode );
    M_Nodes->setunknownNodeMap( theNode->id(), theNode );
879
880
}

881
882
void ZeroDimensionalCircuitData::createKnownNode( const Int & id,
                                                  const zeroDimensionalElementVoltageSourcePtr_Type & theElement )
883
{
884
    zeroDimensionalNodeKnownPtr_Type theNode( new ZeroDimensionalNodeKnown( theElement ) );
885
    theNode->setId( id );
886
887
888
    M_Nodes->setnodeList( theNode );
    M_Nodes->setknownNodeList( theNode );
    M_Nodes->setknownNodeMap( theNode->id(), theNode );
889
890
}

891
void ZeroDimensionalCircuitData::createKnownNode( const Int & id )
892
{
893
    zeroDimensionalNodeKnownPtr_Type theNode( new ZeroDimensionalNodeKnown() );
894
    theNode->setId( id );
895
896
897
    M_Nodes->setnodeList( theNode );
    M_Nodes->setknownNodeList( theNode );
    M_Nodes->setknownNodeMap( theNode->id(), theNode );
898
}
899
900
void ZeroDimensionalCircuitData::fixBC( bcPtr_Type bc )
{
901
902
    M_bc = bc;

903
904
905
906
907
    // Set BC handler in source elements

    ptrVecZeroDimensionalElementCurrentSourcePtr_Type currentElementList = M_Elements->currentSourceList();
    for ( iterZeroDimensionalElementCurrentSource_Type theElement = currentElementList->begin(); theElement != currentElementList->end(); theElement++ )
    {
908
        ( *theElement )->setBC( bc );
909
    }
910
911
912
    ptrVecZeroDimensionalElementVoltageSourcePtr_Type voltageElementList = M_Elements->voltageSourceList();
    for ( iterZeroDimensionalElementVoltageSourcePtr_Type theElement = voltageElementList->begin(); theElement != voltageElementList->end(); theElement++ )
    {
913
        ( *theElement )->setBC( bc );
914
915
916
    }
}

917
918
919
920
void ZeroDimensionalCircuitData::updateCircuitDataFromY( const Real& t,
                                                         const Epetra_Vector* x,
                                                         const Epetra_Vector* x_dot )
{
921
922
923
    const ptrVecZeroDimensionalNodeUnknownPtr_Type& unknownNodeList = M_Nodes->unknownNodeList();

    //update unknown Nodes from solution.
924
925
926
    for ( iterZeroDimensionalNodeUnknown_Type theNode = unknownNodeList ->begin(); theNode != unknownNodeList->end(); theNode++ )
    {
        const Int& variableIndex = ( *theNode )->variableIndex();
927
928
        ( *theNode )->setVoltage( ( *x )[variableIndex] );
        ( *theNode )->setDeltaVoltage( ( *x_dot )[variableIndex] );
929
930
931
932
    }

    //unpdate inductor unknown current
    const ptrVecZeroDimensionalElementPassiveInductorPtr_Type& inductorList = M_Elements->inductorList();
933
934
935
    for ( iterZeroDimensionalElementPassiveInductor_Type theInductor = inductorList ->begin(); theInductor != inductorList->end(); theInductor++ )
    {
        const Int& variableIndex = ( *theInductor )->variableIndex();
936
937
        ( *theInductor )->setCurrent( ( *x )[variableIndex] );
        ( *theInductor )->setDeltaCurrent( ( *x_dot )[variableIndex] );
938
939
940
941
    }

    //update BCs by time
    const ptrVecZeroDimensionalNodeKnownPtr_Type& knownNodeList = M_Nodes->knownNodeList();
942
943
    for ( iterZeroDimensionalNodeKnown_Type theNode = knownNodeList ->begin(); theNode != knownNodeList->end(); theNode++ )
    {
944
945
        ( *theNode )->setVoltageByTime( t );
        ( *theNode )->setDeltaVoltageByTime( t );
946
947
    }
    const ptrVecZeroDimensionalElementCurrentSourcePtr_Type& currentElementList = M_Elements->currentSourceList();
948
949
    for ( iterZeroDimensionalElementCurrentSource_Type theElement = currentElementList ->begin(); theElement != currentElementList->end(); theElement++ )
    {
950
        ( *theElement )->setCurrentByTime( t );
951
952
    }
}
953
void ZeroDimensionalCircuitData::extractSolutionFromY( const Real& t,
954
955
                                                  const Epetra_Vector& y,
                                                  const Epetra_Vector& yp )
956
{
957
958
    //First invoke the updateCircuitDataFromY method (light update)
    updateCircuitDataFromY( t, &y, &yp );
959
960
961

    //deep update (mainly currents), iterate over all elements
    const ptrVecZeroDimensionalElementPtr_Type& elementList = M_Elements->elementList();
962
963
    for ( iterZeroDimensionalElement_Type theElement = elementList ->begin(); theElement != elementList->end(); theElement++ )
    {
964
        ( *theElement )->extractSolution( *M_Nodes );
965
966
967
968
    }

    // Now we can compute voltage source current
    const ptrVecZeroDimensionalElementVoltageSourcePtr_Type& voltageList = M_Elements->voltageSourceList();
969
970
971
972
    for ( iterZeroDimensionalElementVoltageSourcePtr_Type theElement = voltageList ->begin(); theElement != voltageList->end(); theElement++ )
    {
        ( *theElement )->calculateCurrent( *M_Nodes,
                                           *M_Elements );
973
974
975
976
    }

    // check the conservation of current at each node.
    const ptrVecZeroDimensionalNodePtr_Type& nodeList = M_Nodes->nodeList();
977
978
979
    for ( iterZeroDimensionalNode_Type theNode = nodeList ->begin(); theNode != nodeList->end(); theNode++ )
    {
        ( *theNode )->calculateCurrentBalance( *M_Elements );
980
981
982
    }

}
983
984
985
986
void ZeroDimensionalCircuitData::updateABC( matrix_Type& A,
                                            matrix_Type& B,
                                            vector_Type& C )
{
987
988
    A.openCrsMatrix();
    B.openCrsMatrix();
989
990
991