bondgraph-graphviz
来源：互联网 发布：淘宝虚拟禁售商品列表 编辑：程序博客网 时间：2024/06/01 16:31
//******************************************************************************************************
//* gcc -Wall `pkg-config libgvc --cflags --libs` g.c -lgvc -lcgraph
//******************************************************************************************************
/*
* BondGraph.cpp
* Copyright 2009 Jean-Francois Dupuis.
*
* This file is part of libBondGraph.
*
* libBondGraph is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* libBondGraph 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libBondGraph. If not, see <http://www.gnu.org/licenses/>.
*
* This file was created by Jean-Francois Dupuis on 05/03/09.
*/
#include "BondGraph.h"
#include <stdexcept>
#include <cfloat>
#include <stdio.h>
#include <stdlib.h>
#include <fstream>
#include <assert.h>
#include "VectorUtil.h"
#include "BGException.h"
#include "stringutil.h"
using namespace BG;
using namespace std;
BondGraph::BondGraph() : mElementIdCounter(1), mBondIdCounter(1),
mS("S"),mSd("Sd"),mL("L"),mJ("J"),mJj("Jj"),mJIE("Jie"),mJFI("Jfi"),mJFE("Jfe"),mJII("Jii"),
mA("A"), mB("B"), mB2("B2"),mC("C"), mD("D"), mD2("D2"),mCr("Cr"), mDr("Dr"), mD2r("D2r"),
mDynamics(this, &mSimulationLog), mAllowDifferentialCausality(true)
{
//Initialize cound variable
mUcount = -1;
mScount = -1;
mSDcount = -1;
mRcount = -1;
mJcount = -1;
}
BondGraph::~BondGraph()
{
for(unsigned int i = 0; i < mBonds.size(); ++i) {
delete mBonds[i];
mBonds[i] = 0;
}
for(unsigned int i = 0; i < mComponents.size(); ++i) {
delete mComponents[i];
mComponents[i] = 0;
}
}
BondGraph::BondGraph(const BondGraph& inBondGraph) : mDynamics(this, &mSimulationLog) {
*this = inBondGraph;
}
BondGraph& BondGraph::operator=(const BondGraph& inBondGraph) {
//throw runtime_error("BondGraph::operator= not defined!");
//dirty way by lazyness
std::ostringstream lStream;
PACC::XML::Streamer lStreamer(lStream);
inBondGraph.write(lStreamer);
std::istringstream lStreami(lStream.str());
PACC::XML::Document lXMLParser;
lXMLParser.parse(lStreami);
PACC::XML::ConstIterator lRootNode = lXMLParser.getFirstRoot();
this->read(lRootNode);
return *this;
}
/*! \brief Add a component to the bond graph.
* This method add a Component to the elements list. It will attribute the ID
* of the passed component based on the current value of the mElementIdCounter
* and will increment that counter. It will resize the vector if the element
* ID is greater than the size of the vector.
* \param inComponent Component to be added.
*/
void BondGraph::addComponent(Component* inComponent) {
if(inComponent->getId() < 0) {
inComponent->setId(mElementIdCounter++);
}
if( mComponents.size() < inComponent->getId() ) {
mComponents.resize(inComponent->getId());
if(mElementIdCounter <= mComponents.size())
mElementIdCounter = mComponents.size()+1;
}
mComponents[inComponent->getId()-1] = inComponent;
}
/*! \brief Add a bond to the bond graph.
* This method add a Bond to the bond list. It will attribute the ID
* of the passed bond based on the current value of the mBondIdCounter
* and will then increment that counter. It will resize the vector if the bond
* ID is greater than the size of the vector.
* \param inBond Bond to be added
*/
void BondGraph::addBond(Bond* inBond) {
if(inBond->getId() < 0) {
inBond->setId(mBondIdCounter++);
}
if( mBonds.size() < inBond->getId() ) {
mBonds.resize(inBond->getId());
if(mBondIdCounter <= mBonds.size())
mBondIdCounter = mBonds.size()+1;
}
mBonds[inBond->getId()-1] = inBond;
}
/*! \brief Connect two components.
* This method define the connectivity between two Component. The power direction
* is from inOrigin to inDest. The causality is the one defined at the destination component.
* Therefore, if the causality is said to be eEffortCausal, the causality stroke will be applied
* to the inDest side of the bond. Using this method, a new bond is created.
* \param inOrigin Component at the origin of the bond
* \param inDest Component at the end of the bond
* \param inCausality Causality relation between the component, looking at inDest causality.
* \param inBondId Id of the bond that is going to be created.
* \return The newly created bond.
*/
Bond* BondGraph::connect(Component *inOrigin, Component *inDest, Causality inCausality, int inBondId) {
//Create a new bond
Bond *lBond = new Bond(inBondId);
addBond(lBond);
//Create both port
Port *lFromPort = new Port(false);
Port *lToPort = new Port(true);
//Assign causality
if(inCausality == eFlowCausal) {
lFromPort->setCausality(eEffortCausal);
lToPort->setCausality(eFlowCausal);
} else if(inCausality == eEffortCausal) {
lFromPort->setCausality(eFlowCausal);
lToPort->setCausality(eEffortCausal);
}
//Connect the elements together
lBond->connect(lFromPort,lToPort);
inOrigin->connect(lFromPort);
inDest->connect(lToPort);
return lBond;
}
/*! \brief Remove any causality assignment.
* This method set the causality of every ports of all the bond graph components to acausal.
*/
void BondGraph::clearCausality() {
for(unsigned int i = 0; i < mComponents.size(); ++i) {
if(mComponents[i]!=0)
mComponents[i]->clearCausality();
}
}
/*! \brief Propagate causality assignment through the bond graph.
* This recursive methods propagate the causality assignment according to
* the causality constrains of the component. The graph is walked using
* the Port::mVisited flag. When this method assign a causality to a port, it
* will set is flag to true. Hence, the walk throught the graph will not continue
* to a visited port. The walk will stop when causality
* can't be assigned any further according to the component contrains.
* \param inComponent Current component
* \param inPort Incoming port to the current component
*/
void BondGraph::propagateCausality(Component *inComponent, Port* inPort) {
//Apply causality constrain
inComponent->applyCausalityConstrain(inPort);
//Go through all the port of the component
for(vector<Port*>::const_iterator lIter = inComponent->getPorts().begin(); lIter != inComponent->getPorts().end(); ++lIter) {
if((*lIter)->wasVisited())
continue;
//Get the bond connected to this port
Bond* lBond = (*lIter)->getBond();
//Get the other end port
Port* lToPort = lBond->getOtherEnd(*lIter);
//Get the other component
Component* lComponent = lToPort->getComponent();
//Assign both side of the bond the respecting causality
if( (*lIter)->getCausality() == eEffortCausal ) {
if(lToPort->getCausality() != eAcausal && lToPort->getCausality() != eFlowCausal)
throw InvalidCausalityException("Conflicting causality assignment");
lToPort->setCausality(eFlowCausal);
} else if( (*lIter)->getCausality() == eFlowCausal ){
if(lToPort->getCausality() != eAcausal && lToPort->getCausality() != eEffortCausal)
throw InvalidCausalityException("Conflicting causality assignment");
lToPort->setCausality(eEffortCausal);
} else {
return;
}
//Mark visited port
(*lIter)->setVisited();
lToPort->setVisited();
//Go to the next component if the component is not a terminal
if(lComponent->getPorts().size() > 1)
propagateCausality(lComponent, lToPort);
}
}
/*! \brief Assign causality relation to all bonds.
* This methods will compute the causality relation using the sequential
* causal assignment procedure which is describe in System Dynamic, Modeling
* and Simulation of Mechatronic Systems by Karnopp, Margolis and Rosenberg.
*/
void BondGraph::assignCausality() {
//Reset causality and visited port flag
clearCausality();
resetPortFlag();
//Set the causality on the sources
vector<Source*> lSources;
for(unsigned int i = 0; i < mComponents.size(); ++i) {
if(mComponents[i] != 0) {
if(mComponents[i]->getElementType() == BondGraphElement::eSource) {
lSources.push_back( dynamic_cast<Source *>(mComponents[i]) );
}
}
}
//Assign causality on each source and extend the causal implications through the graph. Assignation and propagation
//is done in two step in order to detect conflict with previously assigned source.
for(unsigned int i = 0; i < lSources.size(); ++i) {
lSources[i]->setCausality(); //Set the causality of the source according to his type.
}
for(unsigned int i = 0; i < lSources.size(); ++i) {
propagateCausality(lSources[i], lSources[i]->getPorts()[0]);
}
//Set the causality on C and I
for(unsigned int i = 0; i < mComponents.size(); ++i) {
if(mComponents[i] != 0) {
if(mComponents[i]->getElementType() == BondGraphElement::ePassive) {
Passive* lPassive = dynamic_cast<Passive*> (mComponents[i]);
if(lPassive->getType() == Passive::eCapacitor ||
lPassive->getType() == Passive::eCapacitorM) {
if((lPassive->getPorts())[0]->getCausality()==eAcausal) {
lPassive->getPorts()[0]->setCausality(eFlowCausal);
propagateCausality(lPassive, lPassive->getPorts()[0]);
}
}
else if(lPassive->getType() == Passive::eInductor ||
lPassive->getType() == Passive::eInductorM) {
if(lPassive->getPorts()[0]->getCausality()==eAcausal) {
lPassive->getPorts()[0]->setCausality(eEffortCausal);
propagateCausality(lPassive, lPassive->getPorts()[0]);
}
}
}
}
}
//Assign causality arbitrary on the unassigned elements
for(unsigned int i = 0; i < mComponents.size(); ++i) {
if(mComponents[i] != 0) {
for(unsigned int j = 0; j < mComponents[i]->getPorts().size(); ++j) {
if( mComponents[i]->getPorts()[j]->getCausality()==eAcausal ) {
mComponents[i]->getPorts()[j]->setCausality(eEffortCausal);
propagateCausality(mComponents[i], mComponents[i]->getPorts()[j]);
}
}
}
}
}
/*! \brief Reset port visited flag
* This method marks all Port as unvisited. This is called before to start
* walking the bond graph.
*/
void BondGraph::resetPortFlag() {
for(unsigned int i = 0; i < mComponents.size(); ++i) {
if(mComponents[i] != 0) {
for(vector<Port*>::iterator lIter = mComponents[i]->getPorts().begin(); lIter != mComponents[i]->getPorts().end(); ++lIter) {
(*lIter)->setVisited(false);
}
}
}
}
#ifndef WITHOUT_GRAPHVIZ
/*! \brief Plot a bond graph connection.
* This recursive methods go through all the bond graph connection and add
* an edge to the graph passed as argument. The graph is walked using
* the Port::mVisited flag.
* \param inComponent Current component
* \param inGraph Graphviz graph to write to.
* \param inNodes List of graphic nodes.
* \param inShowBondId Write bond ID on the graph.
*/
void BondGraph::plotGraphConnection(Component* inComponent, Agraph_t *inGraph, std::vector<Agnode_t*> &inNodes, bool inShowBondId) {
for(vector<Port*>::const_iterator lIter = inComponent->getPorts().begin(); lIter != inComponent->getPorts().end(); ++lIter) {
if((*lIter)->wasVisited())
continue;
//Get the bond connected to this port
Bond* lBond = (*lIter)->getBond();
//Get the other end port
Port* lToPort = lBond->getOtherEnd(*lIter);
//Get the other component
Component* lComponent = lToPort->getComponent();
//Add the link to the graph
Agedge_t *lEdge = agedge(inGraph, inNodes[inComponent->getId()-1], inNodes[lComponent->getId()-1],0,0);
if(inShowBondId) {
//Add the bond id as the graph edge label
stringstream lNameStr;
lNameStr << lBond->getId();
char* lName = new char[lNameStr.str().size()+1];
strcpy(lName,lNameStr.str().c_str());
agsafeset(lEdge, "label", lName, "");
delete [] lName;
}
//Check the power direction and causality
if(lToPort->isPowerSide()) {
switch(lToPort->getCausality()) {
case eEffortCausal:
agsafeset(lEdge, "arrowhead", "normal", "");
break;
case eFlowCausal:
agsafeset(lEdge, "arrowhead", "halfopen", "");
break;
default:
agsafeset(lEdge, "arrowhead", "halfopen", "");
//agsafeset(lEdge, "arrowhead", "vee", "");
break;
}
switch((*lIter)->getCausality()) {
case eEffortCausal:
agsafeset(lEdge, "arrowtail", "tee", "");
break;
case eFlowCausal:
agsafeset(lEdge, "arrowtail", "none", "");
break;
default:
agsafeset(lEdge, "arrowtail", "none", "");
//agsafeset(lEdge, "arrowtail", "odot", "");
break;
}
} else {
switch(lToPort->getCausality()) {
case eEffortCausal:
agsafeset(lEdge, "arrowhead", "tee", "");
break;
case eFlowCausal:
agsafeset(lEdge, "arrowhead", "none", "");
break;
default:
agsafeset(lEdge, "arrowhead", "none", "");
//agsafeset(lEdge, "arrowtail", "odot", "");
break;
}
switch((*lIter)->getCausality()) {
case eEffortCausal:
agsafeset(lEdge, "arrowtail", "normal", "");
break;
case eFlowCausal:
agsafeset(lEdge, "arrowtail", "halfopen", "");
break;
default:
agsafeset(lEdge, "arrowtail", "halfopen", "");
//agsafeset(lEdge, "arrowtail", "vee", "");
break;
}
}
//Mark visited port
(*lIter)->setVisited();
lToPort->setVisited();
//Go to the next component if the component is not a terminal
if(lComponent->getPorts().size() > 1)
plotGraphConnection(lComponent, inGraph, inNodes, inShowBondId);
}
}
Component* BondGraph::findUnvisitedComponent() {
for(unsigned int i = 0; i < mComponents.size(); ++i) {
if(mComponents[i] != 0) {
for(vector<Port*>::const_iterator lIter = mComponents[i]->getPorts().begin(); lIter != mComponents[i]->getPorts().end(); ++lIter) {
if((*lIter)->wasVisited())
continue;
return mComponents[i];
}
}
}
return 0;
}
/*! \brief Plot the bond graph.
* This method create a graphical reprensentation of the bond graph using the Graphviz
* library. The default format is svg, but all graphviz supported format can be pass to
* the function.
* \param inFilename The filename where to save the svg bong graph representation
* \param inShowBondId Write bond ID on the graph.
* \param inFormat Output format, default svg
*/
void BondGraph::plotGraph(std::string inFilename, bool inShowBondId, char* inFormat) {
if(mComponents.size() <= 0) {
std::cout << "No elements in graph, cannot be plot" << std::endl;
return;
}
GVC_t* lContext = gvContext();
// Agraph_t *lGraph = agopen("lGraph", AGDIGRAPH,0);
Agraph_t *lGraph = agopen("lGraph", Agdirected,0);
agsafeset(lGraph, "rankdir", "LR", "");
//Create the nodes
vector<Agnode_t*> lNodes(mComponents.size());
for(unsigned int i = 0; i < lNodes.size(); ++i) {
if(mComponents[i] != 0) {
const string& lNameStr = mComponents[i]->getName();
char* lName = new char[lNameStr.size()+1];
strcpy(lName,lNameStr.c_str());
lNodes[i] = agnode(lGraph, lName,0);
agsafeset(lNodes[i], "shape", "none", "");
const string& lDisplayedNameStr = mComponents[i]->getDisplayedName();
char* lDisplayedName = new char[lDisplayedNameStr.size()+1];
strcpy(lDisplayedName,lDisplayedNameStr.c_str());
agsafeset(lNodes[i], "label",lDisplayedName, "");
delete [] lName;
}
}
//Recursively pass through the graph
resetPortFlag();
Component* lComponent;;
while( lComponent = findUnvisitedComponent() ) {
plotGraphConnection(lComponent, lGraph, lNodes, inShowBondId);
}
//char *lFilename = mktemp("libBondGraph.svg.XXXXXX");
char* lName = new char[inFilename.size()+1];
strcpy(lName,inFilename.c_str());
gvLayout(lContext, lGraph, "neato");
gvRenderFilename(lContext, lGraph, inFormat, lName);
//std::cout << "Plot bond graph to " << lName << std::endl;
gvFreeLayout(lContext, lGraph);
agclose(lGraph);
gvFreeContext(lContext);
delete [] lName;
}
#endif
void BondGraph::removeBond(Bond* inBond) {
inBond->getFromPort()->getComponent()->disconnect(inBond->getFromPort());
inBond->getToPort()->getComponent()->disconnect(inBond->getToPort());
//Delete bond
vector<Bond*>::iterator lBondIter = find(mBonds.begin(),mBonds.end(),inBond);
assert(lBondIter != mBonds.end());
*lBondIter = 0;
delete inBond;
}
void BondGraph::removeComponent(std::vector<Component*>& inComponents) {
for(unsigned int i = 0; i < inComponents.size(); ++i) {
removeComponent(inComponents[i]);
}
}
void BondGraph::removeComponent(Component* inComponent) {
vector<Component*>::iterator lComponentIter = find(mComponents.begin(), mComponents.end(), inComponent);
assert(lComponentIter != mComponents.end());
for(unsigned int i = 0; i < inComponent->getPorts().size(); ++i) {
Port* lPort = inComponent->getPorts()[i];
Bond* lBond = lPort->getBond();
//Disconnect the attached component
Port* lOldPort = lBond->getOtherEnd(lPort);
lOldPort->getComponent()->disconnect(lOldPort);
//Check if the bond is part of the output
vector<Bond*>::iterator lBondIter = find(mEffortOutputBonds.begin(),mEffortOutputBonds.end(),lBond);
if( lBondIter != mEffortOutputBonds.end() ) {
*lBondIter = 0;
mEffortOutputBonds.erase(lBondIter);
}
lBondIter = find(mFlowOutputBonds.begin(),mFlowOutputBonds.end(),lBond);
if( lBondIter != mFlowOutputBonds.end() ) {
*lBondIter = 0;
mFlowOutputBonds.erase(lBondIter);
}
//Delete bond
lBondIter = find(mBonds.begin(),mBonds.end(),lBond);
assert(lBondIter != mBonds.end());
*lBondIter = 0;
delete lBond;
lBond = 0;
}
//Once all the bond are deleted, delete the component
delete inComponent;
*lComponentIter = 0;
inComponent = 0;
}
void BondGraph::reconnectComponent(Port* inComponentPort, Component* inNewOrigin) {
Bond* lBond = inComponentPort->getBond();
//Remove connection to the old origin component
Port* lMovingPort = lBond->getOtherEnd(inComponentPort);
lMovingPort->getComponent()->disconnect(lMovingPort);
//Add a new connection to the new origin component
inNewOrigin->connect(lMovingPort);
}
void BondGraph::replaceComponent(BG::Component*& inOldComponent, BG::Component*& inNewComponent) {
vector<Component*>::iterator lComponentIter = find(mComponents.begin(), mComponents.end(), inOldComponent);
assert(lComponentIter != mComponents.end());
//Connect the new component to the same ports
std::vector<Port*> lPorts = inOldComponent->getPorts();
for(unsigned int i = 0; i < lPorts.size(); ++i) {
inNewComponent->connect(lPorts[i]);
}
int lOldId = inOldComponent->getId();
delete inOldComponent;
*lComponentIter = 0;
inOldComponent = 0;
//Add the component
inNewComponent->setId(lOldId);
addComponent(inNewComponent);
}
/*! \brief Simpify Bond Graph
* Remove the redundant componant.
*/
void BondGraph::simplify() {
do{
do{
simplifyRemoveTerminalJunction();
simplifyRemoveEmptyJunction();
simplifyRemoveRedundantJunction();
simplifyRemoveUnconnectedComponent();
} while(simplifyRemoveDoubleBond());
} while( simplifyRemoveRedundantComponents() );
}
void BondGraph::simplifyRemoveUnconnectedComponent() {
//Check if the component is not completely disconnected from the rest. If it's the case, then delete the component.
for(unsigned int i = 0; i < mComponents.size(); ++i) {
Component* lComponent = mComponents[i];
if(lComponent != 0) {
if( lComponent->getPorts().size() == 0 ) {
removeComponent(lComponent);
}
}
}
}
void BondGraph::simplifyRemoveTerminalJunction() {
//Simplify empty junction
for(unsigned int i = 0; i < mComponents.size(); ++i) {
Component* lComponent = mComponents[i];
if(lComponent != 0) {
if( lComponent->getElementType() == BondGraphElement::eJunction ) {
Junction* lJunction = dynamic_cast<Junction*>(lComponent);
//Remove junction if the number of connected component is 2
if(lJunction->getPorts().size() == 1) {
removeComponent(lJunction);
simplifyRemoveTerminalJunction();
return;
}
}
}
}
}
void BondGraph::simplifyRemoveEmptyJunction() {
//Find all junctions
vector<Junction*> lJunctions;
for(unsigned int i = 0; i < mComponents.size(); ++i) {
Component* lComponent = mComponents[i];
if(lComponent != 0) {
if( lComponent->getElementType() == BondGraphElement::eJunction )
lJunctions.push_back(dynamic_cast<Junction*>(lComponent));
}
}
int lRemainingJunction = lJunctions.size();
if(lJunctions.size() > 1) {
for(unsigned int i = 0; i < lJunctions.size() && lRemainingJunction > 1; ++i) {
Junction* lJunction = lJunctions[i];
//Remove junction if the number of connected component is 2
if(lJunction->getPorts().size() == 2) {
//Check that the power is flowing in the same direction
Port* lPort1 = lJunction->getPorts()[0];
Port* lPort2 = lJunction->getPorts()[1];
Component* lComponent1 = lPort1->getBond()->getOtherEnd(lPort1)->getComponent();
Component* lComponent2 = lPort2->getBond()->getOtherEnd(lPort2)->getComponent();
//Check if at least one of the two components is a jonction
if(lComponent1->getElementType() == BondGraphElement::eJunction || lComponent2->getElementType() == BondGraphElement::eJunction) {
if(lPort1->isPowerSide() != lPort2->isPowerSide()) {
//Get the bond connected to the first port
Bond* lBond1 = lPort1->getBond();
//Get the other end port
Port* lToPort1 = lBond1->getOtherEnd(lPort1);
//Get the bond connected to the first port
Bond* lBond2 = lPort2->getBond();
//Get the other end port
Port* lToPort2 = lBond2->getOtherEnd(lPort2);
//Get the component
//Component* lComponent2 = lToPort2->getComponent();
//By pass the second bond
lComponent2->disconnect(lToPort2);
Port* lNewPort = new Port(lToPort2->isPowerSide(), lToPort2->getCausality());
//lNewPort->connect(lBond1);
//lNewPort->connect(lComponent2);
lComponent2->connect(lNewPort);
if(lToPort2->isPowerSide())
lBond1->connect(lToPort1, lNewPort);
else
lBond1->connect(lNewPort, lToPort1);
//Check if the bond is part of the output
vector<Bond*>::iterator lBondIter = find(mEffortOutputBonds.begin(),mEffortOutputBonds.end(),lBond2);
if( lBondIter != mEffortOutputBonds.end() ) {
*lBondIter = lBond1;
}
lBondIter = find(mFlowOutputBonds.begin(),mFlowOutputBonds.end(),lBond2);
if( lBondIter != mFlowOutputBonds.end() ) {
*lBondIter = lBond1;
}
//Delete the junction and bond2
lBondIter = find(mBonds.begin(),mBonds.end(),lBond2);
assert(lBondIter != mBonds.end());
*lBondIter = 0;
delete lBond2;
vector<Component*>::iterator lComponentIter = find(mComponents.begin(), mComponents.end(), lJunction);
assert(lComponentIter != mComponents.end());
*lComponentIter = 0;
delete lJunction;
--lRemainingJunction;
}
}
}
}
}
}
/*! \brief Remove redundant bond between two component
* The method remove the double bond between component. Those generally occur
* after simplifying an empty loop. If the direction are the same, the extra one
* is simply removed. If the direction are not the same, both bond are removed.
* This results in a cut of the bond graph
*/
bool BondGraph::simplifyRemoveDoubleBond() {
bool lFoundDoubleBond = false;
for(unsigned int i = 0; i < mBonds.size(); ++i) {
if(mBonds[i] == 0)
continue;
const Component* lFromComponent = mBonds[i]->getFromPort()->getComponent();
const Component* lToComponent = mBonds[i]->getToPort()->getComponent();
if( lFromComponent == lToComponent ) {
removeBond(mBonds[i]);
lFoundDoubleBond = true;
continue;
}
for(unsigned int j = i+1; j < mBonds.size(); ++j) {
if(mBonds[j] == 0)
continue;
Component* lFromComponent2 = mBonds[j]->getFromPort()->getComponent();
Component* lToComponent2 = mBonds[j]->getToPort()->getComponent();
if( (lFromComponent == lFromComponent2) && (lToComponent == lToComponent2) ) {
//Double bond, same origin, same destination
removeBond(mBonds[i]);
lFoundDoubleBond = true;
j = mBonds.size(); //Continue with the outer for loop
} else if( (lFromComponent == lToComponent2) && (lToComponent == lFromComponent2) ) {
//Link between two components, but with opposite direction.
removeBond(mBonds[i]);
removeBond(mBonds[j]);
lFoundDoubleBond = true;
j = mBonds.size(); //Continue with the outer for loop
}
}
}
return lFoundDoubleBond;
}
void BondGraph::simplifyRemoveRedundantJunction() {
for(unsigned int i = 0; i < mComponents.size(); ++i) {
Component* lComponent = mComponents[i];
if(lComponent != 0) {
//Find all the junction
if( lComponent->getElementType() == BondGraphElement::eJunction ) {
Junction* lJunction = dynamic_cast<Junction*>(lComponent);
Junction::JunctionType lType = lJunction->getType();
for(unsigned int j = 0; j < lJunction->getPorts().size(); ++j) {
//Get the componant connected to each bonds connected to this junction
Bond* lBond = lJunction->getPorts()[j]->getBond();
Port* lIncomingPort = lBond->getOtherEnd(lJunction->getPorts()[j]);
Component* lComponent2 = lIncomingPort->getComponent();
if(lComponent2 == lComponent)
continue;
//Check if the other end of the bond is a junction
if(lComponent2->getElementType() == BondGraphElement::eJunction) {
Junction* lJunction2 = dynamic_cast<Junction*>(lComponent2);
//Check if it's the same type
if(lJunction2->getType() == lType) {
//There is 2 junction of the same type connected together
//their connected components should be all connected to the same junction
vector<Port*> lPortsToReconnect = lJunction2->getPorts();
for(unsigned int k = 0; k < lPortsToReconnect.size(); ++k) {
if(lPortsToReconnect[k] != lIncomingPort) {
Port* lPort = lPortsToReconnect[k];
reconnectComponent(lPort->getBond()->getOtherEnd(lPort), lJunction);
}
}
//remove the extra junction once all te component are reconnected
removeComponent(lJunction2);
}
}
}
}
}
}
}
/*! \brief Simply Redundant Components
* Simplify the bond graph by combining passive component located at the same junction.
* \return True if component are grouped together so that junction could be further simplified
*/
bool BondGraph::simplifyRemoveRedundantComponents() {
bool lSingleComponentIdentified = false;
//Simplify redundant component at the same junction
for(unsigned int i = 0; i < mComponents.size(); ++i) {
if(mComponents[i] != 0) {
Component* lJctComponent = mComponents[i];
if( lJctComponent->getElementType() == BondGraphElement::eJunction ) {
Junction* lJunction = dynamic_cast<Junction*>(lJctComponent);
vector<Passive*> lCapacitors;
vector<Passive*> lResistors;
vector<Passive*> lInductors;
//Find all R,C and I
for(unsigned int j = 0; j < lJunction->getPorts().size(); ++j) {
Bond* lBond = lJunction->getPorts()[j]->getBond();
Component* lComponent = lBond->getOtherEnd(lJunction->getPorts()[j])->getComponent();
if(lComponent->getElementType() == BondGraphElement::ePassive) {
//Exclude the one that are part of the output
if( find(mEffortOutputBonds.begin(),mEffortOutputBonds.end(),lBond) == mEffortOutputBonds.end()
&& find(mFlowOutputBonds.begin(),mFlowOutputBonds.end(),lBond) == mFlowOutputBonds.end() ) {
Passive* lPassive = dynamic_cast<Passive*>(lComponent);
switch(lPassive->getType()) {
case Passive::eResistor:
lResistors.push_back(lPassive);
break;
case Passive::eCapacitor:
lCapacitors.push_back(lPassive);
break;
case Passive::eInductor:
lInductors.push_back(lPassive);
break;
default:
break;
}
}
}
}
//Apply simplification rule and replace the components by a single one
if(lResistors.size() > 1) {
Passive* lNewResistor = combineRedundantPassive(lResistors, lJunction->getType());
addComponent(lNewResistor);
for(unsigned int k = 0; k < lResistors.size(); ++k) {
removeComponent(lResistors[k]);
}
if(lNewResistor != 0)
connect(lJunction, lNewResistor);
if(lJunction->getPorts().size() <= 2)
lSingleComponentIdentified = true;
}
if(lCapacitors.size() > 1) {
Passive* lNewCapacitor = combineRedundantPassive(lCapacitors, lJunction->getType());
addComponent(lNewCapacitor);
for(unsigned int k = 0; k < lCapacitors.size(); ++k) {
removeComponent(lCapacitors[k]);
}
if(lNewCapacitor != 0)
connect(lJunction, lNewCapacitor);
if(lJunction->getPorts().size() <= 2)
lSingleComponentIdentified = true;
}
if(lInductors.size() > 1) {
Passive* lNewInductor = combineRedundantPassive(lInductors, lJunction->getType());
addComponent(lNewInductor);
for(unsigned int k = 0; k < lInductors.size(); ++k) {
removeComponent(lInductors[k]);
}
if(lNewInductor != 0)
connect(lJunction, lNewInductor);
if(lJunction->getPorts().size() <= 2)
lSingleComponentIdentified = true;
}
}
}
}
return lSingleComponentIdentified;
}
Passive* BondGraph::combineRedundantPassive(std::vector<Passive*>& inComponents, Junction::JunctionType inJunctionType) {
if(inComponents.size() == 0)
return 0;
if(inComponents.size() == 1)
return inComponents[0];
Passive* lNewPassive = new Passive(inComponents[0]->getType());
double lValue;
if(inJunctionType == Junction::eZero) {
//Parrallel network simplification
if(lNewPassive->getType() == Passive::eCapacitor) {
lValue = 0;
for(unsigned int i = 0; i < inComponents.size(); ++i) {
if(inComponents[i]->getPorts()[0]->isPowerSide())
lValue += inComponents[i]->getValue();
else
lValue -= inComponents[i]->getValue();
}
} else {
lValue = 1/inComponents[0]->getValue();
for(unsigned int i = 1; i < inComponents.size(); ++i) {
if(inComponents[i]->getPorts()[0]->isPowerSide())
lValue += 1/inComponents[i]->getValue();
else
lValue -= 1/inComponents[i]->getValue();
}
lValue = 1/lValue;
}
} else {
//Serie network simplification
if(lNewPassive->getType() == Passive::eCapacitor) {
lValue = 1/inComponents[0]->getValue();
for(unsigned int i = 1; i < inComponents.size(); ++i) {
if(inComponents[i]->getPorts()[0]->isPowerSide())
lValue += 1/inComponents[i]->getValue();
else
lValue -= 1/inComponents[i]->getValue();
}
lValue = 1/lValue;
} else {
lValue = 0;
for(unsigned int i = 0; i < inComponents.size(); ++i) {
if(inComponents[i]->getPorts()[0]->isPowerSide())
lValue += inComponents[i]->getValue();
else
lValue -= inComponents[i]->getValue();
}
}
}
lNewPassive->setValue(lValue);
return lNewPassive;
}
//Source* BondGraph::combineRedundantSource(vector<Source*>& inComponents, Junction::JunctionType inJunctionType) {
// if(inComponents.size() == 0)
// return 0;
// if(inComponents.size() == 1)
// return inComponents[0];
//
// Source* lNewSource = new Source(inComponents[0]->getType());
// double lValue = 0;
// if(inJunctionType == Junction::eZero) {
// if(lNewSource->getType() != Source::eFlow)
// throw BGException("Multiple effort sources found on a zero junction.");
//
// for(unsigned int i = 0; i < inComponents.size(); ++i) {
// if(inComponents[i]->getPorts()[0]->isPowerSide())
// lValue += inComponents[i]->getValue();
// else
// lValue -= inComponents[i]->getValue();
// }
// } else {
// if(lNewSource->getType() != Source::eEffort)
// throw BGException("Multiple flow sources found on an one junction.");
//
// for(unsigned int i = 0; i < inComponents.size(); ++i) {
// if(inComponents[i]->getPorts()[0]->isPowerSide())
// lValue += inComponents[i]->getValue();
// else
// lValue -= inComponents[i]->getValue();
// }
// }
//
// lNewSource->setValue(lValue);
// return lNewSource;
//}
string BondGraph::serialize(bool inIndent) const {
std::ostringstream lOSS;
PACC::XML::Streamer lStreamer(lOSS, 4);
write(lStreamer, inIndent);
return lOSS.str();
}
void BondGraph::write(PACC::XML::Streamer& inStream, bool inIndent) const {
inStream.openTag("BondGraph", inIndent);
//Write output variable
ostringstream lOSS;
if(mEffortOutputBonds.size() > 0) {
if(mEffortOutputBonds[0] != 0)
lOSS << mEffortOutputBonds[0]->getId();
for(unsigned int i = 1; i < mEffortOutputBonds.size(); ++i) {
if(mEffortOutputBonds[i] != 0)
lOSS << "," << mEffortOutputBonds[i]->getId();
}
string lEffortOutputId = lOSS.str();
inStream.insertAttribute("eout", lEffortOutputId);
}
lOSS.str("");
if(mFlowOutputBonds.size() > 0) {
if(mFlowOutputBonds[0] != 0)
lOSS << mFlowOutputBonds[0]->getId();
for(unsigned int i = 1; i < mFlowOutputBonds.size(); ++i) {
if(mFlowOutputBonds[i] != 0)
lOSS << "," << mFlowOutputBonds[i]->getId();
}
string lFlowOutputId = lOSS.str();
inStream.insertAttribute("fout", lFlowOutputId);
}
//Write component
for(vector<Component*>::const_iterator lComponent = mComponents.begin(); lComponent != mComponents.end(); ++lComponent) {
if(*lComponent != 0)
(*lComponent)->write(inStream,inIndent);
}
for(vector<Bond*>::const_iterator lBond = mBonds.begin(); lBond != mBonds.end(); ++lBond) {
if(*lBond != 0 )
(*lBond)->write(inStream,inIndent);
}
inStream.closeTag();
}
void BondGraph::read(PACC::XML::ConstIterator inNode)
{
//Remove all previous components and bonds
for(unsigned int i = 0; i < mBonds.size(); ++i) {
if(mBonds[i] != 0) {
delete mBonds[i];
mBonds[i] = 0;
}
}
for(unsigned int i = 0; i < mComponents.size(); ++i) {
if(mComponents[i] != 0) {
delete mComponents[i];
mComponents[i] = 0;
}
}
mComponents.resize(0);
mBonds.resize(0);
mEffortOutputBonds.resize(0);
mFlowOutputBonds.resize(0);
if((inNode->getType() != PACC::XML::eData) || (inNode->getValue() != "BondGraph"))
throw BGException("tag <BondGraph> expected!");
//Get output Id
string lEffortOut = inNode->getAttribute("eout");
string lFlowOut = inNode->getAttribute("fout");
//Go through all component
for(PACC::XML::ConstIterator lChild=inNode->getFirstChild(); lChild; lChild=lChild->getNextSibling()) {
if(lChild->getValue() == "Component") {
int lSubType = atoi(lChild->getAttribute("subtype").c_str());
int lId = atoi(lChild->getAttribute("id").c_str());
if(lChild->getAttribute("type") == "Junction") {
addComponent(new Junction(Junction::JunctionType(lSubType),lId));
}
else if(lChild->getAttribute("type") == "Source"){
addComponent(new Source(Source::SourceType(lSubType),atof(lChild->getAttribute("value").c_str()),lId));
}
else if(lChild->getAttribute("type") == "Passive"){
addComponent(new Passive(Passive::PassiveType(lSubType),atof(lChild->getAttribute("value").c_str()),lId));
} else {
throw BGException("Unknown component type!");
}
}
}
mElementIdCounter = mComponents.size()+1;
//Add the bond
for(PACC::XML::ConstIterator lChild=inNode->getFirstChild(); lChild; lChild=lChild->getNextSibling()) {
if(lChild->getValue() == "Bond") {
int lFrom = atoi(lChild->getAttribute("from").c_str());
int lTo = atoi(lChild->getAttribute("to").c_str());
Causality lCausality = Causality(atoi(lChild->getAttribute("causality").c_str()));
int lId = atoi(lChild->getAttribute("id").c_str());
connect(mComponents[lFrom-1],mComponents[lTo-1], lCausality, lId);
}
}
//Adjust bond counter
mBondIdCounter = mBonds.size()+1;
//Assign output bonds
std::istringstream lISS(lEffortOut);
if(!lEffortOut.empty()) {
vector<unsigned int> lEffortOutputId;
lISS >> lEffortOutputId;
for(unsigned int i = 0; i < lEffortOutputId.size(); ++i) {
assert(mBonds[lEffortOutputId[i]-1]->getId() == lEffortOutputId[i]);
//Bond* lBond = mBonds[lEffortOutputId[i]-1];
mEffortOutputBonds.push_back(mBonds[lEffortOutputId[i]-1]);
}
}
if(!lFlowOut.empty()) {
lISS.clear();
lISS.str(lFlowOut);
vector<unsigned int> lFlowOutputId;
lISS >> lFlowOutputId;
for(unsigned int i = 0; i < lFlowOutputId.size(); ++i) {
assert(mBonds[lFlowOutputId[i]-1]->getId() == lFlowOutputId[i]);
mFlowOutputBonds.push_back(mBonds[lFlowOutputId[i]-1]);
}
}
}
/*! \brief Compute the state equation of the bond graph.
* This method will compute the state equation of the bond graph in the form
* \f$\dot{x} = Ax + Bu\f$. If the bond graph contain differential causality
* of the storage component, the form \f$\dot{x} = Ax + Bu + B_2\dot{u}\f$ will
* be used.
*/
void BondGraph::computeStateEquation() {
initializeMatrix();
computeSMatrix();
computeLMatrix();
computeJMatrix();
computeABMatrix();
computeCDMatrix();
computeReducedOutputMatrix();
}
void BondGraph::countComponentType() {
//Count the number of element in each group
mUcount = 0; // Number of source bond
mScount = 0; // Number of storage bond with integral causality
mSDcount = 0; // Number of storage bond with derivative causality
mRcount = 0; // Number of dissipative bond
mJcount = 0; // Number of junction bond
for(vector<Bond*>::iterator lIter = mBonds.begin(); lIter != mBonds.end(); ++lIter) {
if(*lIter != 0) {
(*lIter)->findGroup();
switch( (*lIter)->getGroup() ) {
case Bond::eU:
++mUcount;
break;
case Bond::eS:
++mScount;
break;
case Bond::eSd:
++mSDcount;
break;
case Bond::eR:
++mRcount;
break;
case Bond::eJ:
++mJcount;
break;
default:
throw BGException("BondGraph::countComponentType() : Unknown group type");
}
}
}
}
/*! \brief Initialize the state equations matrices.
* This methods count the number of bond in each Bond::BondGroup. It will resize
* the state equation related matrices according to the size of each group. Each
* bond will be assign an index and an offset according to their position in their
* respective vector.
*
* The \f$x\f$ vector is composed of the state variable corresponding to the
* storage component in integral causality. That is, the charge and momentum variable
* for capacitor and inertia respectively. \f$z\f$ is the coenergy variable vector. It
* form of the effort variable for capacitor and flow variable for inertia. \f$x_d\f$
* and \f$z_d\f$ are the same as previously, but for component in differential
* causality. \f$d_i\f$ and \f$d_o\f$ are the input and output to the dissipative component respectively.
* \f$u\f$ and \f$v\f$ are the output and input of the source component respectively.
*
* Then, the matrices are form such that :
* \f[x = S z\f]
* \f[z_d = S_d x_d\f]
* \f[d_o = L d_i\f]
* \f[ \left[\begin{array}{c}\dot{x}\\z_d\\d_i\\v\end{array}\right] = J \left[\begin{array}{c}z\\\dot{x_d}\\d_o\\u\end{array}\right]\f]
* \f[ \left[\begin{array}{c} \dot{x} \\ z_d \\ d_i \\ v \end{array}\right] = J_{FE} \left[\begin{array}{c} z \\ \dot{x_d} \\ d_o \\ u \end{array}\right] \f]
* \f[ \left[\begin{array}{c} J_{e} \\ J_{f} \end{array}\right] = J_{IE} \left[\begin{array}{c} z \\ \dot{x_d} \\ d_o \\ u \end{array}\right] \f]
* \f[ \left[\begin{array}{c} \dot{x} \\ z_d \\ d_i \\ v \end{array}\right] = J_{FI} \left[\begin{array}{c} J_{e} \\ J_{f} \end{array}\right] \f]
* \f[ \left[\begin{array}{c} J_{e} \\ J_{f} \end{array}\right] = J_{II} \left[\begin{array}{c} J_{e} \\ J_{f} \end{array}\right] \f]
* where \f$J_e\f$ and \f$J_f\f$ are the effort and flow vector of the junctions bonds. For example \f$[J_{e},J_{f}]^T = [e_1,e_2,e_3,f_1,f_2,f_3]^T\f$
*
* The index of the bond give the position of the bond into is respective vector and the
* the offset give the additional index into the concatenated vector.
*/
void BondGraph::initializeMatrix() {
//Count the number of element in each group
countComponentType();
int lUSRcount = mUcount + mScount + mSDcount + mRcount;
//Resize the matrix according to the relevant group size
mS.setZero(mScount,mScount);
mSd.setZero(mSDcount,mSDcount);
mL.setZero(mRcount,mRcount);
mJ.setZero(lUSRcount,lUSRcount);
mJFE.setZero(lUSRcount,lUSRcount);
mJIE.setZero(2*mJcount,lUSRcount);
mJFI.setZero(lUSRcount,2*mJcount);
mJII.setZero(2*mJcount,2*mJcount);
//Assign matrix index to each bond
int lXidx = 0;
int lXDidx = 0;
int lDidx = 0;
int lUidx = 0;
int lJidx = 0;
int lSidx = 0;
for(vector<Bond*>::iterator lIter = mBonds.begin(); lIter != mBonds.end(); ++lIter) {
if(*lIter != 0) {
switch( (*lIter)->getGroup() ) {
case Bond::eS:
if((*lIter)->isEffortInput()) {
(*lIter)->setMatrixIndex(lXidx, 0, lXidx, -lXidx-1);
} else {
(*lIter)->setMatrixIndex(lXidx, 0, -lXidx-1, lXidx);
}
++lXidx;
(*lIter)->setStorageIndex(lSidx++);
break;
case Bond::eSd:
if((*lIter)->isEffortInput()) {
(*lIter)->setMatrixIndex(lXDidx, mScount, lXDidx+mScount+mSDcount, lXDidx+mScount);
} else {
(*lIter)->setMatrixIndex(lXDidx, mScount, lXDidx+mScount, lXDidx+mScount+mSDcount);
}
++lXDidx;
(*lIter)->setStorageIndex(lSidx++);
break;
case Bond::eR:
if((*lIter)->isEffortInput()) {
(*lIter)->setMatrixIndex(lDidx, mScount+mSDcount, lDidx+mScount+mSDcount+mSDcount+mRcount, lDidx+mScount+mSDcount+mSDcount);
} else {
(*lIter)->setMatrixIndex(lDidx, mScount+mSDcount, lDidx+mScount+mSDcount+mSDcount, lDidx+mScount+mSDcount+mSDcount+mRcount);
}
++lDidx;
break;
case Bond::eU:
if((*lIter)->isEffortInput()) {
(*lIter)->setMatrixIndex(lUidx, mScount+mSDcount+mRcount, -lUidx-mScount-1, lUidx+mScount+mSDcount+mSDcount+mRcount+mRcount);
} else {
(*lIter)->setMatrixIndex(lUidx, mScount+mSDcount+mRcount, lUidx+mScount+mSDcount+mSDcount+mRcount+mRcount, -lUidx-mScount-1);
}
++lUidx;
break;
case Bond::eJ:
(*lIter)->setMatrixIndex(lJidx, 0, lJidx+mScount+mSDcount+mSDcount+mRcount+mRcount+mUcount, lJidx+mScount+mSDcount+mSDcount+mRcount+mRcount+mUcount+mJcount, mJcount);
++lJidx;
break;
default:
throw BGException("BondGraph::initializeMatrix() : Unknown bond group");
}
#ifdef DEBUG_EQN
if( (*lIter)->getGroup() == Bond::eJ ) {
cout << "Bond " << (*lIter)->getId() << " is part of group " << (*lIter)->getGroup() << " with index " << (*lIter)->getMatrixIndex().index << " and joffset " << (*lIter)->getMatrixIndex().joffset << ". Effort and flow index are : " << (*lIter)->getMatrixIndex().effort << ", " << (*lIter)->getMatrixIndex().flow << endl;
} else {
cout << "Bond " << (*lIter)->getId() << " is part of group " << (*lIter)->getGroup() << " with index " << (*lIter)->getMatrixIndex().index << " and offset " << (*lIter)->getMatrixIndex().index << ". Effort and flow index are : " << (*lIter)->getMatrixIndex().effort << ", " << (*lIter)->getMatrixIndex().flow << endl;
}
#endif
}
}
//Compute the U vector and name vector
//mU.resize(mUcount);
mSources.resize(mUcount);
mNameU.resize(mUcount);
mNameX.resize(mScount);
mNameXd.resize(mSDcount);
for(unsigned int i = 0; i < mComponents.size(); ++i) {
if(mComponents[i] != 0) {
switch( mComponents[i]->getComponentGroup() ) {
case Component::eU:
mSources[mComponents[i]->getPorts()[0]->getBond()->getMatrixIndex().index] = dynamic_cast<Source*>(mComponents[i]);//->getValue();
mNameU[mComponents[i]->getPorts()[0]->getBond()->getMatrixIndex().index] = mComponents[i]->getName();
break;
case Component::eS:
if(mComponents[i]->getPorts()[0]->getBond()->getGroup() != Bond::eSd) {
mNameX[mComponents[i]->getPorts()[0]->getBond()->getMatrixIndex().index] = mComponents[i]->getName();
} else {
mNameXd[mComponents[i]->getPorts()[0]->getBond()->getMatrixIndex().index] = mComponents[i]->getName();
}
break;
default:
break;
}
}
}
//Form the output name vector
mOutputName.resize(2*mBonds.size() - mUcount);
for(vector<Bond*>::iterator lIter = mBonds.begin(); lIter != mBonds.end(); ++lIter) {
if(*lIter != 0) {
Bond::MatrixIndex lIdx = (*lIter)->getMatrixIndex();
if(lIdx.effort >= 0)
mOutputName[lIdx.effort] = string((*lIter)->getName()+".e");
if(lIdx.flow >= 0)
mOutputName[lIdx.flow] = string((*lIter)->getName()+".f");
}
}
// std::cout << "NameX: " << mNameX << std::endl;
// std::cout << "NameXd: " << mNameXd << std::endl;
// std::cout << "mOutputName: " << mOutputName << std::endl;
}
/*! \brief Compute the S matrix
* This method compute the diagonal S matrix describes in #initializeMatrix.
*/
void BondGraph::computeSMatrix() {
for(vector<Component*>::iterator lIter = mComponents.begin(); lIter != mComponents.end(); ++lIter) {
if(*lIter != 0) {
if( (*lIter)->getComponentGroup() == Component::eS) {
if( (*lIter)->getPorts()[0]->getBond()->getGroup() == Bond::eS ) {
int lIdx = (*lIter)->getPorts()[0]->getBond()->getMatrixIndex().index;
if( (*lIter)->getValue() == 0 ) {
mS(lIdx,lIdx) = DBL_MAX;
} else {
mS(lIdx,lIdx) = 1/(*lIter)->getValue();
}
} else {
int lIdx = (*lIter)->getPorts()[0]->getBond()->getMatrixIndex().index;
mSd(lIdx,lIdx) = (*lIter)->getValue();
}
}
}
}
#ifdef DEBUG_EQN
PACC::XML::Streamer lStream(cout);
mS.write(lStream); cout << endl;
mSd.write(lStream); cout << endl;
#endif
}
/*! \brief Compute the L matrix
* This method compute the diagonal L matrix describes in #initializeMatrix.
*/
void BondGraph::computeLMatrix() {
for(vector<Component*>::iterator lIter = mComponents.begin(); lIter != mComponents.end(); ++lIter) {
if(*lIter != 0) {
if( (*lIter)->getComponentGroup() == Component::eR) {
int lIdx = (*lIter)->getPorts()[0]->getBond()->getMatrixIndex().index;
if( (*lIter)->getPorts()[0]->getCausality() == eEffortCausal ) {
if( (*lIter)->getValue() == 0 ) {
mL(lIdx,lIdx) = DBL_MAX;
} else {
mL(lIdx,lIdx) = 1/(*lIter)->getValue();
}
} else {
mL(lIdx,lIdx) = (*lIter)->getValue();
}
}
}
}
#ifdef DEBUG_EQN
PACC::XML::Streamer lStream(cout);
mL.write(lStream); cout << endl;
#endif
}
/*! \brief Compute the J matrix
* First this methods fill the elements of the JIE,JFI,JFE and JII matrix from the bond graph junction
* structure. Then, the J matrix is computed as \f$J = J_{FE} + J_{FI}(I-J_{II})^{-1} J_{IE} \f$
* if \f$J_{II}\f$ is not zero or \f$J = J_{FE} + J_{FI}J_{IE} \f$ otherwise.
*/
void BondGraph::computeJMatrix() {
//Compute JE,JF,JFE and JI matrix
for(vector<Component*>::iterator lIter = mComponents.begin(); lIter != mComponents.end(); ++lIter) {
if(*lIter != 0) {
if( (*lIter)->getComponentGroup() == Component::eJ ) {
if( (*lIter)->getElementType() == BondGraphElement::eJunction ) {
Junction* lJunction = dynamic_cast<Junction*> (*lIter);
if( lJunction->getType() == Junction::eZero ) {
//Find the base, i.e. the effort causal connected port
Port* lBasePort = 0;
for(vector<Port*>::iterator lIter = lJunction->getPorts().begin(); lIter != lJunction->getPorts().end(); ++lIter) {
if( (*lIter)->getCausality() == eEffortCausal ) {
lBasePort = *lIter;
break;
}
}
if(lBasePort == 0) {
throw BGException(string("Invalid causality assignment. No baseport found on junction ")+lJunction->getName());
}
//Fill the matrix according to the equation ei = eb and the sum relation for f
for(vector<Port*>::iterator lIter = lJunction->getPorts().begin(); lIter != lJunction->getPorts().end(); ++lIter) {
if( (*lIter)->getCausality() != eEffortCausal ) {
setJMatrixElement(lBasePort, *lIter, lJunction->getType());
}
}
} else { //Junction eOne
//Find the base, i.e. the flow causal connected port
Port* lBasePort = 0;
for(vector<Port*>::iterator lIter = lJunction->getPorts().begin(); lIter != lJunction->getPorts().end(); ++lIter) {
if( (*lIter)->getCausality() == eFlowCausal ) {
lBasePort = *lIter;
break;
}
}
if(lBasePort == 0) {
throw BGException(string("Invalid causality assignment. No baseport found on junction ")+lJunction->getName());
}
//Fill the matrix according to the equation fi = fb and the sum relation for e
for(vector<Port*>::iterator lIter = lJunction->getPorts().begin(); lIter != lJunction->getPorts().end(); ++lIter) {
if( (*lIter)->getCausality() != eFlowCausal ) {
setJMatrixElement(lBasePort, *lIter, lJunction->getType());
}
}
}
} else { //GY or TR
Passive* lPassive = dynamic_cast<Passive*> (*lIter);
setJMatrixElement(lPassive->getPorts(),lPassive->getType());
}
}
}
}
#ifdef DEBUG_EQN
PACC::XML::Streamer lStream(cout);
mJFE.write(lStream); cout << endl;
mJIE.write(lStream); cout << endl;
mJFI.write(lStream); cout << endl;
mJII.write(lStream); cout << endl;
#endif
//Compute the J matrix with J = JFE + JFI*inv(I-JI)*JIE if JII is not zero or J = JFE + JFI*JIE if JII is zero.
bool lIsJIZero = true;
for(unsigned i = 0; i < mJII.getRows(); ++i) {
for(unsigned int j = 0; j < mJII.getCols(); ++j) {
if(mJII(i,j) != 0) {
lIsJIZero = false;
break;
}
}
}
if(!lIsJIZero) {
PACC::Matrix lEye;
lEye.setIdentity(2*mJcount);
PACC::Matrix lInvJI = (lEye-mJII).invert();
mJj = lInvJI*mJIE;
mJ = mJFE + mJFI*lInvJI*mJIE;
} else {
mJj = mJIE;
mJ = mJFE + mJFI*mJIE;
}
#ifdef DEBUG_EQN
mJ.write(lStream); cout << endl;
#endif
}
//#ifdef DEBUG_EQN
//void BondGraph::displayEquation(const std::vector< Port* > &inPorts, Passive::PassiveType inType) {
// if(inType == Passive::eTransformer || inType == Passive::eTransformerM) {
// cout << "Equation : e" << inPorts[1]->getBond()->getId() << " = e" << inPorts[0]->getBond()->getId() << endl;
// } else {
//
// }
//}
//#endif
/*! \brief Assign element of the appropriate J matrix according to the constrains of the TF and GY elements
* This methods fill the JFE, JIE, JFI or JII according to the type of bond involved.
* \param inPorts Ports of the component
* \param inType Type of component
*/
void BondGraph::setJMatrixElement(const std::vector< Port* > &inPorts, Passive::PassiveType inType) {
assert(inPorts.size() == 2); //There should have to port to TF or GY
Port *lBasePort, *lPort;
//Find the base port
if(inPorts[0]->isPowerSide()) {
assert(!inPorts[1]->isPowerSide()); //Power should flow
lBasePort = inPorts[0];
lPort = inPorts[1];
} else {
assert(inPorts[1]->isPowerSide()); //Power should flow
lBasePort = inPorts[1];
lPort = inPorts[0];
}
//Find the value of the transformation
bool lIsEffortCausal;
double lValue;
if( lBasePort->getCausality() == eEffortCausal ) {
lValue = 1/lBasePort->getComponent()->getValue();
lIsEffortCausal = true;
}
else {
lValue = lBasePort->getComponent()->getValue();
lIsEffortCausal = false;
}
//If the power direction is the same -> -1, if opposed -> 1
if( (lBasePort->isPowerSide() && lPort->isPowerSide()) ||
(!lBasePort->isPowerSide() && !lPort->isPowerSide()) ) {
lValue = -lValue;
}
int lIdx1 = lBasePort->getBond()->getMatrixIndex().index + lBasePort->getBond()->getMatrixIndex().offset;
int lIdx2 = lPort->getBond()->getMatrixIndex().index + lPort->getBond()->getMatrixIndex().offset;
if( lBasePort->getBond()->getGroup() == Bond::eS || lBasePort->getBond()->getGroup() == Bond::eSd || lBasePort->getBond()->getGroup() == Bond::eR || lBasePort->getBond()->getGroup() == Bond::eU) {
if( lPort->getBond()->getGroup() == Bond::eS || lPort->getBond()->getGroup() == Bond::eSd || lPort->getBond()->getGroup() == Bond::eR || lPort->getBond()->getGroup() == Bond::eU) {
mJFE(lIdx1, lIdx2) = lValue;
mJFE(lIdx2, lIdx1) = lValue;
} else { //SUR and J
int lIdx2F = lIdx2 + lPort->getBond()->getMatrixIndex().joffset;
if(inType == Passive::eTransformer || inType == Passive::eTransformerM) {
if(lIsEffortCausal) {
mJFI(lIdx1,lIdx2F) = lValue;
mJIE(lIdx2,lIdx1) = lValue;
} else {
mJFI(lIdx1,lIdx2) = lValue;
mJIE(lIdx2F,lIdx1) = lValue;
}
} else {
if(lIsEffortCausal) {
mJFI(lIdx1, lIdx2) = lValue;
mJIE(lIdx2F, lIdx1) = lValue;
} else {
mJFI(lIdx1, lIdx2F) = lValue;
mJIE(lIdx2, lIdx1) = lValue;
}
}
}
} else {
if( lPort->getBond()->getGroup() == Bond::eS || lPort->getBond()->getGroup() == Bond::eSd || lPort->getBond()->getGroup() == Bond::eR || lPort->getBond()->getGroup() == Bond::eU) {
//J and USR
int lIdx1F = lIdx1 + lBasePort->getBond()->getMatrixIndex().joffset;
if(inType == Passive::eTransformer || inType == Passive::eTransformerM) {
if(lIsEffortCausal) {
mJIE(lIdx1F,lIdx2) = lValue;
mJFI(lIdx2,lIdx1) = lValue;
} else {
mJIE(lIdx1,lIdx2) = lValue;
mJFI(lIdx2,lIdx1F) = lValue;
}
} else {
if(lIsEffortCausal) {
mJIE(lIdx1F, lIdx2) = lValue;
mJFI(lIdx2, lIdx1) = lValue;
} else {
mJIE(lIdx1, lIdx2) = lValue;
mJFI(lIdx2, lIdx1F) = lValue;
}
}
} else {//J and J
int lIdx1F = lIdx1 + lBasePort->getBond()->getMatrixIndex().joffset;
int lIdx2F = lIdx2 + lPort->getBond()->getMatrixIndex().joffset;
if(inType == Passive::eTransformer || inType == Passive::eTransformerM) {
if(lIsEffortCausal) {
mJII(lIdx1F,lIdx2F) = lValue;
mJII(lIdx2,lIdx1) = lValue;
} else {
mJII(lIdx1,lIdx2) = lValue;
mJII(lIdx2F,lIdx1F) = lValue;
}
} else {
if(lIsEffortCausal) {
mJII(lIdx1F, lIdx2) = lValue;
mJII(lIdx2F, lIdx1) = lValue;
} else {
mJII(lIdx1, lIdx2F) = lValue;
mJII(lIdx2, lIdx1F) = lValue;
}
}
}
}
}
/*! \brief Assign element of the appropriate J matrix according to the constrains of the junction-0 or junction-1.
* Set the element of the JFE, JIE, JFI or JII matrix based on the type bonds involved.
* \param inBasePort The base bond of the junction
* \param inPort A dependent bond of the junction
* \param inType The type of junction
*/
void BondGraph::setJMatrixElement(const Port* inBasePort, const Port* inPort, Junction::JunctionType inType) {
#ifdef DEBUG_EQN
displayEquation(inBasePort,inPort,inType);
#endif
int lIdxBase = inBasePort->getBond()->getMatrixIndex().index + inBasePort->getBond()->getMatrixIndex().offset;
int lIdx = inPort->getBond()->getMatrixIndex().index + inPort->getBond()->getMatrixIndex().offset;
if( inBasePort->getBond()->getGroup() == Bond::eS || inBasePort->getBond()->getGroup() == Bond::eSd || inBasePort->getBond()->getGroup() == Bond::eR || inBasePort->getBond()->getGroup() == Bond::eU) {
if( inPort->getBond()->getGroup() == Bond::eS || inPort->getBond()->getGroup() == Bond::eSd || inPort->getBond()->getGroup() == Bond::eR || inPort->getBond()->getGroup() == Bond::eU) {
//Set equality relation
mJFE(lIdx, lIdxBase) = 1;
#ifdef DEBUG_EQN
cout << "Jfe(" << lIdx << ", " << lIdxBase << ") = 1" << endl;
#endif
//Set sum equality relation. If the power direction is the same -> -1, if opposed -> 1
if( (inBasePort->isPowerSide() && inPort->isPowerSide()) ||
(!inBasePort->isPowerSide() && !inPort->isPowerSide()) ) {
mJFE(lIdxBase, lIdx) = -1;
#ifdef DEBUG_EQN
cout << "Jfe(" << lIdxBase << ", " << lIdx << ") = -1" << endl;
#endif
} else {
mJFE(lIdxBase, lIdx) = 1;
#ifdef DEBUG_EQN
cout << "Jfe(" << lIdxBase << ", " << lIdx << ") = 1" << endl;
#endif
}
} else {//inPort->getBond()->getGroup() == Bond::eJ
int lIdxE = lIdx;
int lIdxF = lIdx;
if(inType == Junction::eZero) {
lIdxF = lIdx + inPort->getBond()->getMatrixIndex().joffset;
} else {
lIdxE = lIdx + inPort->getBond()->getMatrixIndex().joffset;
}
mJIE(lIdxE, lIdxBase) = 1;
#ifdef DEBUG_EQN
cout << "Je(" << lIdxE << ", " << lIdxBase << ") = 1" << endl;
#endif
//Set sum equality relation. If the power direction is the same -> -1, if opposed -> 1
if( (inBasePort->isPowerSide() && inPort->isPowerSide()) ||
(!inBasePort->isPowerSide() && !inPort->isPowerSide()) ) {
mJFI(lIdxBase, lIdxF) = -1;
#ifdef DEBUG_EQN
cout << "Jf(" << lIdxBase << ", " << lIdxF << ") = -1" << endl;
#endif
} else {
mJFI(lIdxBase, lIdxF) = 1;
#ifdef DEBUG_EQN
cout << "Jf(" << lIdxBase << ", " << lIdxF << ") = 1" << endl;
#endif
}
}
} else { //inBasePort->getBond()->getGroup() == Bond::eJ
if( inPort->getBond()->getGroup() == Bond::eS || inPort->getBond()->getGroup() == Bond::eSd || inPort->getBond()->getGroup() == Bond::eR || inPort->getBond()->getGroup() == Bond::eU) {
int lIdxBaseE = lIdxBase;
int lIdxBaseF = lIdxBase;
if(inType == Junction::eZero) {
lIdxBaseF = lIdxBase + inBasePort->getBond()->getMatrixIndex().joffset; //Change from lIdx
} else {
lIdxBaseE = lIdxBase + inBasePort->getBond()->getMatrixIndex().joffset; //Change from lIdx
}
mJFI(lIdx, lIdxBaseE) = 1;
#ifdef DEBUG_EQN
cout << "Jf(" << lIdx << ", " << lIdxBaseE << ") = 1" << endl;
#endif
//Set sum equality relation. If the power direction is the same -> -1, if opposed -> 1
if( (inBasePort->isPowerSide() && inPort->isPowerSide()) ||
(!inBasePort->isPowerSide() && !inPort->isPowerSide()) ) {
mJIE(lIdxBaseF, lIdx) = -1;
#ifdef DEBUG_EQN
cout << "Je(" << lIdxBaseF << ", " << lIdx << ") = -1" << endl;
#endif
} else {
mJIE(lIdxBaseF, lIdx) = 1;
#ifdef DEBUG_EQN
cout << "Je(" << lIdxBaseF << ", " << lIdx << ") = 1" << endl;
#endif
}
} else {//inPort->->getBond()getGroup() == Bond::eJ
int lIdxBaseE = lIdxBase;
int lIdxBaseF = lIdxBase;
int lIdxE = lIdx;
int lIdxF = lIdx;
if(inType == Junction::eZero) {
lIdxF = lIdx + inPort->getBond()->getMatrixIndex().joffset;
lIdxBaseF = lIdxBase + inPort->getBond()->getMatrixIndex().joffset;//Change from lIdx
} else {
lIdxE = lIdx + inPort->getBond()->getMatrixIndex().joffset;
lIdxBaseE = lIdxBase + inPort->getBond()->getMatrixIndex().joffset;//Change from lIdx
}
mJII(lIdxE, lIdxBaseE) = 1;
#ifdef DEBUG_EQN
cout << "Ji(" << lIdxE << ", " << lIdxBaseE << ") = 1" << endl;
#endif
//Set sum equality relation. If the power direction is the same -> -1, if opposed -> 1
if( (inBasePort->isPowerSide() && inPort->isPowerSide()) ||
(!inBasePort->isPowerSide() && !inPort->isPowerSide()) ) {
mJII(lIdxBaseF, lIdxF) = -1;
#ifdef DEBUG_EQN
cout << "Ji(" << lIdxBaseF << ", " << lIdxF << ") = -1" << endl;
#endif
} else {
mJII(lIdxBaseF, lIdxF) = 1;
#ifdef DEBUG_EQN
cout << "Ji(" << lIdxBaseF << ", " << lIdxF << ") = 1" << endl;
#endif
}
}
}
}
#ifdef DEBUG_EQN
void BondGraph::displayEquation(const Port* inBasePort, const Port* inPort, Junction::JunctionType inType) {
if(inType == Junction::eZero)
cout << "Equation : e" << inPort->getBond()->getId() << " = e" << inBasePort->getBond()->getId() << endl;
else
cout << "Equation : f" << inPort->getBond()->getId() << " = f" << inBasePort->getBond()->getId() << endl;
if( (inBasePort->isPowerSide() && inPort->isPowerSide()) ||
(!inBasePort->isPowerSide() && !inPort->isPowerSide()) ) {
if(inType == Junction::eZero)
cout << "Equation : f" << inBasePort->getBond()->getId() << " = -f" << inPort->getBond()->getId() << " + ... " << endl;
else
cout << "Equation : e" << inBasePort->getBond()->getId() << " = -e" << inPort->getBond()->getId() << " + ... " << endl;
} else {
if(inType == Junction::eZero)
cout << "Equation : f" << inBasePort->getBond()->getId() << " = f" << inPort->getBond()->getId() << " + ... " << endl;
else
cout << "Equation : e" << inBasePort->getBond()->getId() << " = e" << inPort->getBond()->getId() << " + ... " << endl;
}
}
#endif
void BondGraph::computeCDMatrix() {
//Check if there is derivative causality
if(mSDcount == 0) {
//Extract the submatrix from J
PACC::Matrix lJxz("Jxz"), lJxd("Jxd"), lJxu("Jxu");
PACC::Matrix lJdz("Jdz"), lJdd("Jdd"), lJdu("Jdu");
PACC::Matrix lJvz("Jvz"), lJvd("Jvd"), lJvu("Jvu");
if(mScount > 0) {
mJ.extract(lJxz, 0, mScount-1, 0, mScount-1);
if(mRcount > 0)
mJ.extract(lJxd, 0, mScount-1, mScount, mScount+mRcount-1);
if(mUcount > 0)
mJ.extract(lJxu, 0, mScount-1, mScount + mRcount, mScount + mRcount + mUcount-1);
}
if(mRcount > 0) {
mJ.extract(lJdd, mScount, mScount+mRcount-1, mScount, mScount+mRcount-1);
if(mScount > 0)
mJ.extract(lJdz, mScount, mScount+mRcount-1, 0, mScount-1);
if(mUcount > 0)
mJ.extract(lJdu, mScount, mScount+mRcount-1, mScount + mRcount, mScount + mRcount + mUcount-1);
}
if(mUcount > 0) {
mJ.extract(lJvu, mScount+mRcount, mScount + mRcount + mUcount-1, mScount + mRcount, mScount + mRcount + mUcount-1);
if(mScount > 0)
mJ.extract(lJvz, mScount+mRcount, mScount + mRcount + mUcount-1, 0, mScount-1);
if(mRcount > 0)
mJ.extract(lJvd, mScount+mRcount, mScount + mRcount + mUcount-1, mScount, mScount+mRcount-1);
}
//Extract submatrix from Jj
PACC::Matrix lJjz,lJjd,lJju;
if(mJj.size() > 0) {
if(mScount > 0)
mJj.extractColumns(lJjz, 0, mScount-1);
if(mRcount > 0)
mJj.extractColumns(lJjd, mScount, mScount+mRcount-1);
if(mUcount > 0)
mJj.extractColumns(lJju, mScount+mRcount, mScount+mRcount+mUcount-1);
}
//Compute Cdo = L*inv(I-Jdd*L)*Jdz*S
//Compute Cdi = I + Jdd*L*inv(I-Jdd*L)*Jdz*S ///not the same as in the code
//Compute Cj = (Jjz + Jjd*L*inv(I-Jdd*L)*Jdz)*S
//Compute Dj = Jju + Jjd*L*inv(I-Jdd*L)*Jdu
//where Jj = inv(I-Jii)*Jie
PACC::Matrix lCdo,lDdo,lCdi,lDdi,lCj,lDj,lCv,lDv;
mC = mS;
if(mRcount > 0) {
PACC::Matrix lI;
lI.setIdentity(mL.getCols());
PACC::Matrix lInv = mL*((lI - lJdd*mL).invert());
PACC::Matrix lM = (lI + lJdd*lInv);
if(mScount > 0) {
lCdi = lM*lJdz*mS;
lCdo = lInv*lJdz*mS;
mC.concatenateRows(mC,lCdi);
mC.concatenateRows(mC,lCdo);
}
if(mUcount > 0) {
mD.resize(0,0);
mD.resize(mScount,mUcount);
lDdo = lInv*lJdu;
lDdi = lM*lJdu;
if(mScount > 0)
lCv = (lJvz+lJvd*lInv*lJdz)*mS;
lDv = lJvd*lInv*lJdu+lJvu; //To verify : lInv*lJdu+lJvu;
if(mScount > 0)
mC.concatenateRows(mC,lCv);
mD.concatenateRows(mD,lDdi);
mD.concatenateRows(mD,lDdo);
mD.concatenateRows(mD,lDv);
if(mJj.size() > 0) {
lDj = lJju + lJjd*lInv*lJdu;
mD.concatenateRows(mD,lDj);
}
} else {
mD.resize(0,0);
}
if(mJj.size() > 0) {
if(mScount > 0) {
lCj = (lJjz + lJjd*lInv*lJdz)*mS;
mC.concatenateRows(mC,lCj);
}
}
} else {
if(mUcount > 0) {
if(mScount > 0) {
lCv = lJvz*mS;
mC.concatenateRows(mC,lCv);
}
mD.resize(0,0);
mD.resize(mScount,mUcount);
lDv = lJvu;
mD.concatenateRows(mD,lDv); //Problem here when S ->| 0 ->| I
}
if(mJj.size() > 0) {
if(mScount > 0) {
lCj = lJjz*mS;
mC.concatenateRows(mC,lCj);
}
if(mUcount > 0) {
mD.resize(mScount,mUcount);
lDv = lJvu;
mD.concatenateRows(mD,lDv); //to verify
lDj = lJju;
mD.concatenateRows(mD,lDj);
} else {
mD.resize(0,0);
}
}
}
mD2.resize(0,0);
#ifdef DEBUG_EQN
PACC::XML::Streamer lStream(cout);
mC.write(lStream); cout << endl;
mD.write(lStream); cout << endl;
#endif
} else {
//There is some derivative causality
//Extract the submatrix from J
PACC::Matrix Jii("Jii"), Jid("Jid"), Jil("Jil"), Jiu("Jiu");
PACC::Matrix Jdi("Jdi"), Jdu("Jdu");
PACC::Matrix Jli("Jli"), Jld("Jld"), Jll("Jll"), Jlu("Jlu");
PACC::Matrix Jvi("Jvi"), Jvd("Jvd"), Jvl("Jvl"), Jvu("Jvu");
if(mScount > 0) {
mJ.extract(Jii, 0, mScount-1, 0, mScount-1);
if(mSDcount > 0)
mJ.extract(Jdi, mScount, mScount+mSDcount-1, 0, mScount-1);
if(mRcount > 0)
mJ.extract(Jli, mScount+mSDcount, mScount+mSDcount+mRcount-1, 0, mScount-1);
}
if(mSDcount > 0) {
if(mScount > 0)
mJ.extract(Jid, 0, mScount-1, mScount, mScount+mSDcount-1);
if(mRcount > 0)
mJ.extract(Jld, mScount+mSDcount, mScount+mSDcount+mRcount-1, mScount, mScount+mSDcount-1);
}
if(mRcount > 0) {
if(mScount > 0)
mJ.extract(Jil, 0, mScount-1, mScount+mSDcount, mScount+mSDcount+mRcount-1);
mJ.extract(Jll, mScount+mSDcount, mScount+mSDcount+mRcount-1, mScount+mSDcount, mScount+mSDcount+mRcount-1);
}
if(mUcount > 0) {
if(mScount > 0)
mJ.extract(Jiu, 0, mScount-1, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
if(mSDcount > 0)
mJ.extract(Jdu, mScount, mScount+mSDcount-1, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
if(mRcount > 0)
mJ.extract(Jlu, mScount+mSDcount, mScount+mSDcount+mRcount-1, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
}
//Extract submatrix for computing v
if(mUcount > 0) {
unsigned int lIdx1 = mScount+mSDcount+mRcount;
unsigned int lIdx2 = mScount+mSDcount+mRcount+mUcount-1;
mJ.extract(Jvu, lIdx1, lIdx2, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
mJ.extract(Jvd, lIdx1, lIdx2, mScount, mScount+mSDcount-1);
mJ.extract(Jvd, lIdx1, lIdx2, mScount, mScount+mSDcount-1);
if(mScount > 0)
mJ.extract(Jvi, lIdx1, lIdx2, 0, mScount-1);
if(mRcount > 0)
mJ.extract(Jvl, lIdx1, lIdx2, mScount+mSDcount, mScount+mSDcount+mRcount-1);
}
//Extract submatrix from Jj
PACC::Matrix Jji,Jjd,Jju,Jjl;
if(mJj.size() > 0) {
mJj.extractColumns(Jjd, mScount, mScount+mSDcount-1);
if(mScount > 0)
mJj.extractColumns(Jji, 0, mScount-1);
if(mRcount > 0)
mJj.extractColumns(Jjl, mScount+mSDcount, mScount+mSDcount+mRcount-1);
if(mUcount > 0)
mJj.extractColumns(Jju, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
}
//Matrix declaration
PACC::Matrix Czd,Dzd,Dzd2,Cxd,Dxd,Dxd2,Cdo,Ddo,Ddo2,Cdi,Ddi,Ddi2,Cj,Dj,Dj2,Cv,Dv,Dv2;
PACC::Matrix F,G1,G2,H,K,M;
//Useful invert
PACC::Matrix lI; lI.setIdentity(Jll.getRows());
PACC::Matrix P = mL*(lI - Jll).invert();
//Useful substitution
if(mScount != 0) {
F = mSd*Jdi*mS*mA;
G1 = mSd*Jdi*mS*mB;
H = P*(Jli*mS + Jld*F);
G2 = mSd*Jdi*mS*mB2 + mSd*Jdu;
} else {
G1.resize(mSDcount,mUcount);
G2 = mSd*Jdu;
}
K = P*(Jlu + Jld*G1);
M = P*Jld*G2;
//Compute C and D matrix
mC = mS;
if(mScount > 0) {
Czd = Jdi*mS;
mC.concatenateRows(mC,Czd);
Cxd = F;
mC.concatenateRows(mC,Cxd);
}
if(mRcount > 0 && mScount > 0) {
Cdo = H;
Cdi = Jli*mS+Jld*F+Jll*H;
mC.concatenateRows(mC,Cdi);
mC.concatenateRows(mC,Cdo);
}
if(mUcount > 0) {
mD.resize(0,0);
mD.resize(mScount,mUcount);
mD2.resize(0,0);
mD2.resize(mScount,mUcount);
Dzd = Jdu;
Dzd2 = PACC::Matrix(Jdu.getRows(),Jdu.getCols());
mD.concatenateRows(mD,Dzd);
mD2.concatenateRows(mD2,Dzd2);
if(mScount > 0) {
Cv = Jvi*mS+Jvd*F+Jvl*H;
mC.concatenateRows(mC,Cv);
}
Dxd = G1;
Dxd2 = G2;
mD.concatenateRows(mD,Dxd);
mD2.concatenateRows(mD2,Dxd2);
if(mRcount > 0) {
Ddo = K;
Ddo2 = M;
Ddi = Jlu + Jld*G1 + Jll*K;
Ddi2 = Jld*G2 + Jll*M;
mD.concatenateRows(mD,Ddi);
mD2.concatenateRows(mD2,Ddi2);
mD.concatenateRows(mD,Ddo);
mD2.concatenateRows(mD2,Ddo2);
}
Dv = Jvu + Jvd*G1 + Jvl*K;
Dv2 = Jvd*G2 + Jvl*M;
mD.concatenateRows(mD,Dv);
mD2.concatenateRows(mD2,Dv2);
if(mJcount > 0) {
Dj = Jjl*K + Jjd*G1 + Jju;
Dj2 = Jjl*M + Jjd*G2;
mD.concatenateRows(mD,Dj);
mD2.concatenateRows(mD2,Dj2);
}
}
if(mJcount > 0 && mScount > 0) {
Cj = Jji*mS + Jjl*H + Jjd*F;
mC.concatenateRows(mC,Cj);
}
#ifdef DEBUG_EQN
PACC::XML::Streamer lStream(cout);
mC.write(lStream); cout << endl;
mD.write(lStream); cout << endl;
mD2.write(lStream); cout << endl;
cout << getOutputVariableNames() << endl;
#endif
}
//Verification
if(mScount != 0 && mUcount != 0 && (mC.getRows() != mD.getRows()) ) {
throw BGException("Matrix C and D don't have the same number of rows");
}
if(mUcount != 0 && mSDcount != 0 && (mD.getRows() != mD2.getRows()) ) {
throw BGException("Matrix D and D2 don't have the same number of rows");
}
}
/*! \brief Compute the value of the A and B matrix of the state equation
* In the case where there is no differential causality, the state equation
* can be formulated as \f$\dot{x} = Ax + Bu\f$. The matrices are then computed
* as :
* \f[
* \begin{array}{cl}
* A &= (J_{xz} + J_{xd}L(I - J_{dd}L)^{-1}J_{dz})S \\
* B &= J_{xu} + J_{xd}L(I - J_{dd}L)^{-1}J_{du}
* \end{array}
* \f]
* where the submatrices are extracted as :
* \f[ \begin{array}{cl}
* \left[\begin{array}{c}\dot{x}\\d_i\\v\end{array}\right] &= J \left[\begin{array}{c}z\\d_o\\u\end{array}\right] \\
* \left[\begin{array}{c}\dot{x}\\d_i\end{array}\right] &= \left[\begin{array}{ccc} J_{xz} & J_{xd} & J_{xu} \\ J_{dz} & J_{dd} & J_{du} \end{array}\right] \left[\begin{array}{c}z\\\dot{x_d}\\d_o\\u\end{array}\right]
* \end{array}\f]
*
* In the case where there is differential causality involved, the state equation
* are formulated as \f$\dot{x} = Ax + Bu + B_2\dot{u}\f$. The matrices are then
* computed as :
* \f[
* \begin{array}{cl}
* K &= J_{il}L(I - LJ_{ll})^{-1}\\
* T_{ii} &= J_{ii}S + K J_{li}S \\
* T_{iu} &= J_{iu} + K J{lu} \\
* T_f &= I - J_{id} S_{d} J_{di} S - K J_{ld} S_d J_{di} S \\
* T_{i\dot{u}} &= K J_{ld} S_{d} J_{du} + J_{id} S_{d} J_{du} \\
* A &= T_f^{-1}T_{ii} \\
* B &= T_f^{-1}T_{iu} \\
* B_2 &= T_f^{-1}T_{i\dot{u}}
* \end{array}
* \f]
* where the submatrices are extracted as :
* \f[ \begin{array}{cl}
* \left[\begin{array}{c}\dot{x}\\z_d\\d_i\\v\end{array}\right] &= J \left[\begin{array}{c}z\\\dot{x_d}\\d_o\\u\end{array}\right] \\
* \left[\begin{array}{c}\dot{x}\\z_d\\d_i\end{array}\right] &= \left[\begin{array}{cccc} J_{ii} & J_{id} & J_{il} & J_{iu} \\ J_{di} & J_{dd} & J_{dl} & J_{du} \\ J_{li} & J_{ld} & J_{ll} & J_{lu} \end{array}\right] \left[\begin{array}{c}z\\\dot{x_d}\\d_o\\u\end{array}\right]
* \end{array}\f]
*/
void BondGraph::computeABMatrix() {
//Check if there is any state to the system. If not, clear A and B.
if(mScount == 0) {
mA.resize(0,0);
mB.resize(0,0);
#ifdef DEBUG_EQN
PACC::XML::Streamer lStream(cout);
mA.write(lStream); cout << endl;
mB.write(lStream); cout << endl;
#endif
return;
}
//Check if there is derivative causality
if(mSDcount == 0) {
//Extract the submatrix from J
PACC::Matrix lJxz, lJxd, lJxu;
PACC::Matrix lJdz, lJdd, lJdu;
if(mScount > 0)
mJ.extract(lJxz, 0, mScount-1, 0, mScount-1);
if(mRcount > 0)
mJ.extract(lJxd, 0, mScount-1, mScount, mScount+mRcount-1);
if(mUcount > 0)
mJ.extract(lJxu, 0, mScount-1, mScount + mRcount, mScount + mRcount + mUcount-1);
if(mRcount > 0) {
mJ.extract(lJdd, mScount, mScount+mRcount-1, mScount, mScount+mRcount-1);
if(mScount > 0)
mJ.extract(lJdz, mScount, mScount+mRcount-1, 0, mScount-1);
if(mUcount > 0)
mJ.extract(lJdu, mScount, mScount+mRcount-1, mScount + mRcount, mScount + mRcount + mUcount-1);
}
//Compute A and B matrix as A = (Jxz + Jxd*L*inv(I - Jdd*L)*Jdz)*S and B = Jxu + Jxd*L*inv(I - Jdd*L)*Jdu
if(mRcount > 0) {
PACC::Matrix lI;
lI.setIdentity(mL.getCols());
PACC::Matrix lInv = lJxd*mL*(lI - lJdd*mL).invert(); //Jxd*L*inv(I - Jdd*L)
mA = (lJxz + lInv*lJdz)*mS;
if(mUcount > 0)
mB = lJxu + lInv*lJdu;
else
mB.resize(0,0);
}
else {
mA = lJxz*mS;
if(mUcount > 0)
mB = lJxu;
else
mB.resize(0,0);
}
mB2.resize(0,0);
} else {
//Compute using A, B in presence of differential causality
//Extract the submatrix from J
PACC::Matrix Jii("Jii"), Jid("Jid"), Jil("Jil"), Jiu("Jiu");
PACC::Matrix Jdi("Jdi"), Jdu("Jdu");
PACC::Matrix Jli("Jli"), Jld("Jld"), Jll("Jll"), Jlu("Jlu");
if(mScount > 0) {
mJ.extract(Jii, 0, mScount-1, 0, mScount-1);
if(mSDcount > 0)
mJ.extract(Jdi, mScount, mScount+mSDcount-1, 0, mScount-1);
if(mRcount > 0)
mJ.extract(Jli, mScount+mSDcount, mScount+mSDcount+mRcount-1, 0, mScount-1);
}
if(mSDcount > 0) {
mJ.extract(Jid, 0, mScount-1, mScount, mScount+mSDcount-1);
if(mRcount > 0)
mJ.extract(Jld, mScount+mSDcount, mScount+mSDcount+mRcount-1, mScount, mScount+mSDcount-1);
}
if(mRcount > 0) {
mJ.extract(Jil, 0, mScount-1, mScount+mSDcount, mScount+mSDcount+mRcount-1);
mJ.extract(Jll, mScount+mSDcount, mScount+mSDcount+mRcount-1, mScount+mSDcount, mScount+mSDcount+mRcount-1);
}
if(mUcount > 0) {
mJ.extract(Jiu, 0, mScount-1, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
mJ.extract(Jdu, mScount, mScount+mSDcount-1, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
if(mRcount > 0)
mJ.extract(Jlu, mScount+mSDcount, mScount+mSDcount+mRcount-1, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
}
/*
if(mScount > 0) {
mJ.extract(Jii, 0, mScount-1, 0, mScount-1);
mJ.extract(Jdi, mScount, mScount+mSDcount-1, 0, mScount-1);
mJ.extract(Jli, mScount+mSDcount, mScount+mSDcount+mRcount-1, 0, mScount-1);
}
if(mSDcount > 0) {
mJ.extract(Jid, 0, mScount-1, mScount, mScount+mSDcount-1);
mJ.extract(Jld, mScount+mSDcount, mScount+mSDcount+mRcount-1, mScount, mScount+mSDcount-1);
}
if(mRcount > 0) {
mJ.extract(Jil, 0, mScount-1, mScount+mSDcount, mScount+mSDcount+mRcount-1);
mJ.extract(Jll, mScount+mSDcount, mScount+mSDcount+mRcount-1, mScount+mSDcount, mScount+mSDcount+mRcount-1);
}
if(mUcount > 0) {
mJ.extract(Jiu, 0, mScount-1, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
mJ.extract(Jdu, mScount, mScount+mSDcount-1, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
mJ.extract(Jlu, mScount+mSDcount, mScount+mSDcount+mRcount-1, mScount+mSDcount+mRcount, mScount+mSDcount+mRcount+mUcount-1);
}
*/
#ifdef DEBUG_EQN
PACC::XML::Streamer lStream(cout);
Jii.write(lStream); cout << endl;
Jdi.write(lStream); cout << endl;
Jli.write(lStream); cout << endl;
Jid.write(lStream); cout << endl;
Jld.write(lStream); cout << endl;
Jil.write(lStream); cout << endl;
Jll.write(lStream); cout << endl;
Jiu.write(lStream); cout << endl;
Jdu.write(lStream); cout << endl;
Jlu.write(lStream); cout << endl;
#endif
//Check if Jdi = 0 and Jdu = 0, if so Tf.invert = identity and Tiuu = 0
bool lJdiZero = true;
bool lJduZero = true;
for(unsigned int i = 0; i < Jdi.getRows(); ++i) {
for(unsigned int j = 0; j < Jdi.getCols(); ++j) {
if( Jdi(i,j) != 0 ) {
lJdiZero = false;
break;
}
}
}
if(lJdiZero) {
for(unsigned int i = 0; i < Jdu.getRows(); ++i) {
for(unsigned int j = 0; j < Jdu.getCols(); ++j) {
if( Jdu(i,j) != 0 ) {
lJdiZero = false;
break;
}
}
}
}
//Compute A, B and B2 matrices
PACC::Matrix Il;
Il.setIdentity(mRcount);
PACC::Matrix Tii("Tii"), Tiu("Tiu");
PACC::Matrix lInv;
if(mRcount > 0)
lInv = Jil*mL*(Il - mL*Jll).invert();
Tii = Jii*mS + lInv*Jli*mS;
if(mUcount > 0)
Tiu = Jiu + lInv*Jlu;
if( !(lJdiZero && lJduZero) ) {
PACC::Matrix Tf("Tf"), Tiuu("Tiuu"), Ii;
Ii.setIdentity(mScount);
if(mRcount > 0)
Tf = Ii - Jid*mSd*Jdi*mS - lInv*Jld*mSd*Jdi*mS;
else
Tf = Ii - Jid*mSd*Jdi*mS;
PACC::Matrix Tf_inv = Tf.invert();
mA = Tf_inv * Tii;
if(mUcount > 0) {
if(mRcount > 0)
Tiuu = lInv*Jld*mSd*Jdu + Jid*mSd*Jdu;
else
Tiuu = Jid*mSd*Jdu;
//Tiuu = lInv*Jld*mSd*Jdi*mS + Jid*mSd*Jdu;
mB = Tf_inv * Tiu;
mB2 = Tf_inv * Tiuu;
} else {
mB.resize(0,0);
mB2.resize(0,0);
}
#ifdef DEBUG_EQN
Tiuu.write(lStream); cout << endl;
Tf.write(lStream); cout << endl;
#endif
} else {
mA = Tii;
mB = Tiu;
mB2 = PACC::Matrix(mScount,mUcount,0);
}
#ifdef DEBUG_EQN
Tii.write(lStream); cout << endl;
Tiu.write(lStream); cout << endl;
#endif
}
#ifdef DEBUG_EQN
PACC::XML::Streamer lStream(cout);
mA.write(lStream); cout << endl;
mB.write(lStream); cout << endl;
if(mSDcount != 0)
mB2.write(lStream); cout << endl;
#endif
}
void BondGraph::setInitialState(const std::vector<double> &inValues) {
if(mScount != inValues.size()) {
throw BGException("The number of initial values must be equal to the number of states in the system.");
}
mDynamics.setInput(getSourceValues(0));
mDynamics.initialize(inValues);
}
const vector<double>& BondGraph::getInputs() const {
return( mDynamics.getInput() );
}
const vector<double>& BondGraph::getInputsDt() const {
return( mDynamics.getDuDt() );
}
//void BondGraph::getStateVariables(vector<double>& outStates) const {
// outStates = mDynamics.getStateVariables();
//}
const vector<double>& BondGraph::getOutputVariables() {
return mDynamics.getObservableVariables();
}
const vector<double>& BondGraph::getStateVariables() const {
return mDynamics.getStateVariables();
}
//void BondGraph::getOutputVariables(vector<double>& outOutputs) {
// outOutputs = mDynamics.getObservableVariables();
//}
vector<double> BondGraph::getSourceValues(double inTime) const {
vector<double> lOutValues(mSources.size());
for(unsigned int i = 0; i < mSources.size(); ++i) {
lOutValues[i] = mSources[i]->getValue(inTime);
}
return lOutValues;
}
void BondGraph::simulate(double inStopTime, double inSimulationStep/*, double inLogStep*/) {
mDynamics.reset();
mDynamics.setTimeStep(inSimulationStep/10);
while(mDynamics.getTime() < inStopTime) {
//Get source value
mDynamics.setInput( getSourceValues(mDynamics.getTime()) );
//Simulate the system with the source value
mDynamics.simulate(mDynamics.getTime() + inSimulationStep);
mDynamics.writeLog();
}
}
/*! \brief Write the simulation log to file
* \param inFilename Name of the file to write the simulation to.
*/
void BondGraph::writeSimulationLog(std::string inFilename) {
if(mSimulationLog.size() > 0) {
ofstream lStream(inFilename.c_str());
//Write the header
lStream << "#" << mSimulationLog.begin()->first;
for(std::map<string, vector<double> >::const_iterator lIter = ++(mSimulationLog.begin()); lIter != mSimulationLog.end(); ++lIter) {
lStream << "," << lIter->first;
}
lStream << endl;
//Write the data
lStream.precision(15);
for(unsigned int i = 0; i < mSimulationLog.begin()->second.size(); ++i) {
lStream << (mSimulationLog.begin()->second)[i];
for(std::map<string, vector<double> >::const_iterator lIter = ++(mSimulationLog.begin()); lIter != mSimulationLog.end(); ++lIter) {
lStream << "," << (lIter->second)[i];
}
lStream << endl;
}
}
}
/*! \brief Defined the output bonds.
* This methods assign the bonds that will form the output vector.
* \param inEffortBonds Bonds for which the effort is monitored.
* \param inFlowBonds Bonds for which the flow is monitored.
*/
void BondGraph::setOutputBonds(const std::vector<Bond*> &inEffortBonds,const std::vector<Bond*> &inFlowBonds) {
mEffortOutputBonds = inEffortBonds;
mFlowOutputBonds = inFlowBonds;
//computeReducedOutputMatrix();
}
void BondGraph::setOutputBonds(Bond* inEffortBond, Bond* inFlowBond) {
if(inEffortBond != 0)
mEffortOutputBonds.resize(1, inEffortBond);
else
mEffortOutputBonds.resize(0);
if(inFlowBond != 0)
mFlowOutputBonds.resize(1,inFlowBond);
else
mFlowOutputBonds.resize(0);
//computeReducedOutputMatrix();
}
void BondGraph::addOutputBonds(Bond* inBond, ValueType inType) {
if(inType == eEffort) {
vector<Bond*>::iterator lIter = find(mEffortOutputBonds.begin(), mEffortOutputBonds.end(), inBond);
if(lIter == mEffortOutputBonds.end()) {
mEffortOutputBonds.push_back(inBond);
}
} else {
vector<Bond*>::iterator lIter = find(mFlowOutputBonds.begin(), mFlowOutputBonds.end(), inBond);
if(lIter == mFlowOutputBonds.end()) {
mFlowOutputBonds.push_back(inBond);
}
}
}
void BondGraph::computeReducedOutputMatrix() {
if( (mC.size() > 0 || mD.size() > 0) &&
!((mEffortOutputBonds.size() == 0) && (mFlowOutputBonds.size() == 0)) )
{
mOutputName.resize(0);
mCr.resize(0,mC.getCols());
mDr.resize(0,mD.getCols());
mD2r.resize(0,mD2.getCols());
PACC::Matrix lRow;
for(vector<Bond*>::iterator lIter = mEffortOutputBonds.begin(); lIter != mEffortOutputBonds.end(); ++lIter) {
mOutputName.push_back(string((*lIter)->getName()+".e"));
int lIdx = (*lIter)->getMatrixIndex().effort;
if(lIdx >= 0) {
//Variable from the full output vector
if(mC.size() > 0)
mCr.concatenateRows(mCr,mC.extractRow(lRow,lIdx));
if(mD.size() > 0)
mDr.concatenateRows(mDr,mD.extractRow(lRow,lIdx));
} else {
lIdx = -lIdx-1;
if( lIdx < mA.getRows() ) {
//Variable form the dot{x} vector
if(mC.size() > 0)
mCr.concatenateRows(mCr,mA.extractRow(lRow,lIdx));
if(mD.size() > 0)
mDr.concatenateRows(mDr,mB.extractRow(lRow,lIdx));
}
else {
//Variable is part of the u vector
if(mC.size() > 0) {
lRow.resize(0,0);
lRow.resize(1,mA.getCols()); //resize
mCr.concatenateRows(mCr,lRow);
}
if(mD.size() > 0) {
lRow.resize(0,0);
lRow.resize(1,mD.getCols());
lRow(0,lIdx - mA.getRows() ) = 1;
mDr.concatenateRows(mDr,lRow);
};
}
}
if(mC.size() > 0) {
(*lIter)->getMatrixIndex().reffort = mCr.getRows()-1;
} else {
(*lIter)->getMatrixIndex().reffort = mDr.getRows()-1;
}
}
for(vector<Bond*>::iterator lIter = mFlowOutputBonds.begin(); lIter != mFlowOutputBonds.end(); ++lIter) {
mOutputName.push_back(string((*lIter)->getName()+".f"));
int lIdx = (*lIter)->getMatrixIndex().flow;
if(lIdx >= 0) {
//Variable from the full output vector
if(mC.size() > 0)
mCr.concatenateRows(mCr,mC.extractRow(lRow,lIdx));
if(mD.size() > 0)
mDr.concatenateRows(mDr,mD.extractRow(lRow,lIdx));
} else {
//Variable form the dot{x} vector
lIdx = -lIdx-1;
if( lIdx < mA.getRows() ) {
if(mC.size() > 0)
mCr.concatenateRows(mCr,mA.extractRow(lRow,lIdx));
if(mD.size() > 0)
mDr.concatenateRows(mDr,mB.extractRow(lRow,lIdx));
}
else {
//Variable is part of the u vector
if(mC.size() > 0) {
lRow.resize(0,0);
lRow.resize(1,mA.getCols());
mCr.concatenateRows(mCr,lRow);
}
if(mD.size() > 0) {
lRow.resize(0,0);
lRow.resize(1,mD.getCols());
lRow(0,lIdx - mA.getRows()) = 1;
mDr.concatenateRows(mDr,lRow);
}
}
}
if(mC.size() > 0) {
(*lIter)->getMatrixIndex().rflow = mCr.getRows()-1;
} else {
(*lIter)->getMatrixIndex().rflow = mDr.getRows()-1;
}
}
} else {
mCr = mC;
mDr = mD;
mD2r = mD2;
}
#ifdef DEBUG_EQN
PACC::XML::Streamer lStream(cout);
mCr.write(lStream); cout << endl;
mDr.write(lStream); cout << endl;
mD2r.write(lStream); cout << endl;
#endif
}
//Go through all junctions
bool BondGraph::parsingCompare(Component* inComponent1, Component* inComponent2) {
//Compare the number of connection to that component
if( inComponent1->getPorts().size() != inComponent2->getPorts().size() ) return false;
//Compare the type of component connected to the junction
if( !(inComponent1->countConnectionByGroup() == inComponent2->countConnectionByGroup()) ) return false;
//Go through all the port of the component to continue to parse the graph.
for(vector<Port*>::const_iterator lIter1 = inComponent1->getPorts().begin(); lIter1 != inComponent1->getPorts().end(); ++lIter1) {
if((*lIter1)->wasVisited())
continue;
//Get the bond connected to this port
Bond* lBond1 = (*lIter1)->getBond();
//Get the other end port
Port* lToPort1 = lBond1->getOtherEnd(*lIter1);
//Get the other component
Component* lComponent1 = lToPort1->getComponent();
//Mark visited port
(*lIter1)->setVisited();
lToPort1->setVisited();
//Parse only junction
if( lComponent1->getComponentGroup() == Component::eJ ) {
//Here with need to test all branch, if there is no match branch, return false
bool lMatchingBranch = false;
for(vector<Port*>::const_iterator lIter2 = inComponent2->getPorts().begin(); lIter2 != inComponent2->getPorts().end(); ++lIter2) {
//Get the bond connected to this port
Bond* lBond2 = (*lIter2)->getBond();
//Get the other end port
Port* lToPort2 = lBond2->getOtherEnd(*lIter2);
//Get the other component
Component* lComponent2 = lToPort2->getComponent();
if( lComponent2->getComponentGroup() == Component::eJ ) {
if( parsingCompare(lComponent1,lComponent2) ) {
lMatchingBranch = true;
break; //Found a matching branch
}
}
}
if(!lMatchingBranch)
return false;
}
}
return true;
}
/*! \brief Compare the structure of two bond graph
* \param inBondgraph The bond graph to be compared against
* \param inStart1 A bond that serve as a starting point of comparison
* \param inStart2 A bond that serve as a starting point of comparison from the inBondgraph
* \return True is the two bond graph have the same structure
*/
bool BondGraph::compare(BondGraph& inBondgraph, const Bond* inStart1, const Bond* inStart2) {
//The function assume that the bond graph have been already simplified
//Compare the number of each type of elements
countComponentType();
inBondgraph.countComponentType();
if( mUcount != inBondgraph.mUcount) return false;
if( mScount+mSDcount != inBondgraph.mScount+inBondgraph.mSDcount) return false;
if( mRcount != inBondgraph.mRcount) return false;
if( mJcount != inBondgraph.mJcount) return false;
//Clear the visited port flag
resetPortFlag();
if(inStart1 != 0 && inStart2 != 0) {
//Get both starting junction
Component* lToComponent1 = inStart1->getToPort()->getComponent();
Component* lToComponent2 = inStart2->getToPort()->getComponent();
if( lToComponent1->getComponentGroup() != lToComponent2->getComponentGroup() ) return false;
if( lToComponent1->getComponentGroup() == Component::eJ ) {
return parsingCompare(lToComponent1,lToComponent2);
}
else {
//if not a junction, the other one must be one
return parsingCompare(inStart1->getFromPort()->getComponent(),inStart2->getFromPort()->getComponent());
}
return true;
} else {
//There is no starting point, hence all combinaison of starting point should be tried.
//Get the smallest component set to make starting points
std::vector<unsigned int> lCounts;
if(mUcount != 0) lCounts.push_back(mUcount);
if(mScount+mSDcount != 0) lCounts.push_back(mScount+mSDcount);
if(mRcount != 0) lCounts.push_back(mRcount);
if(mJcount != 0) lCounts.push_back(mJcount);
int lMin = *(std::min_element( lCounts.begin(), lCounts.end() ) );
//Go true all the matching pair
for(vector<Component*>::const_iterator lIter1 = mComponents.begin(); lIter1 != mComponents.end(); ++lIter1) {
if( *lIter1 == 0 ) continue;
Component *lComponent1 = *lIter1;
if(mJcount == lMin) {
if( (*lIter1)->getComponentGroup() != Component::eJ )
continue;
} else {
if(mUcount == lMin) { //Use U set
if( (*lIter1)->getComponentGroup() != Component::eU )
continue;
} else if( (mScount+mSDcount) == lMin ) { //Use S set
if( (*lIter1)->getComponentGroup() != Component::eS)
continue;
} else { //Use R set
if( (*lIter1)->getComponentGroup() != Component::eR )
continue;
}
//Get the junction that is connected to the component
lComponent1 = (*lIter1)->getPorts()[0]->getBond()->getOtherEnd((*lIter1)->getPorts()[0])->getComponent();
}
for(vector<Component*>::const_iterator lIter2 = inBondgraph.mComponents.begin(); lIter2 != inBondgraph.mComponents.end(); ++lIter2) {
if( *lIter2 == 0 ) continue;
Component *lComponent2 = *lIter2;
if(mJcount == lMin) {
if( (*lIter2)->getComponentGroup() != Component::eJ )
continue;
} else {
if(mUcount == lMin) { //Use U set
if( (*lIter2)->getComponentGroup() != Component::eU )
continue;
} else if( (mScount+mSDcount) == lMin ) { //Use S set
if( (*lIter2)->getComponentGroup() != Component::eS)
continue;
} else { //Use R set
if( (*lIter2)->getComponentGroup() != Component::eR )
continue;
}
//Get the junction that is connected to the component
lComponent2 = (*lIter2)->getPorts()[0]->getBond()->getOtherEnd((*lIter2)->getPorts()[0])->getComponent();
}
//Compare the pair
if(parsingCompare(lComponent1,lComponent2))
return true;
}
}
return false; //No matching structure found from all the possible combination of the starting set.
}
}
//***************************************************************************************************
//*
//***************************************************************************************************
I'm trying to run a example from the "Using Graphviz as a library" in http://www.graphviz.org/Documentation.php.
#include <gvc.h>int main(int argc, char **argv){ Agraph_t *g; Agnode_t *n, *m; Agedge_t *e; Agsym_t *a; GVC_t *gvc; /* set up a graphviz context */ gvc = gvContext(); /* parse command line args - minimally argv[0] sets layout engine */ gvParseArgs(gvc, argc, argv); /* Create a simple digraph */ g = agopen("g", Agdirected); n = agnode(g, "n", 1); m = agnode(g, "m", 1); e = agedge(g, n, m, 0, 1); /* Set an attribute - in this case one that affects the visible rendering */ agsafeset(n, "color", "red", ""); /* Compute a layout using layout engine from command line args */ gvLayoutJobs(gvc, g); /* Write the graph according to -T and -o options */ gvRenderJobs(gvc, g); /* Free layout data */ gvFreeLayout(gvc, g); /* Free graph structures */ agclose(g); /* close output file, free context, and return number of errors */ return (gvFreeContext(gvc));}I'm compiling and linking with : gcc -Wall pkg-config libgvc --cflags --libs *.c -o EXE -lgvc
and then I see this result:
graph.c: In function ‘main’:graph.c:14:18: error: ‘Agdirected’ undeclared (first use in this function)graph.c:14:18: note: each undeclared identifier is reported only once for each function it appears ingraph.c:15:2: error: too many arguments to function ‘agnode’In file included from /usr/include/graphviz/types.h:717:0, from /usr/include/graphviz/gvc.h:20, from graph.c:1:/usr/include/graphviz/graph.h:185:22: note: declared heregraph.c:16:2: error: too many arguments to function ‘agnode’In file included from /usr/include/graphviz/types.h:717:0, from /usr/include/graphviz/gvc.h:20, from graph.c:1:/usr/include/graphviz/graph.h:185:22: note: declared heregraph.c:17:2: error: too many arguments to function ‘agedge’In file included from /usr/include/graphviz/types.h:717:0, from /usr/include/graphviz/gvc.h:20, from graph.c:1:/usr/include/graphviz/graph.h:192:22: note: declared heregraph.c:7:11: warning: unused variable ‘a’ [-Wunused-variable]graph.c:6:12: warning: variable ‘e’ set but not used [-Wunused-but-set-variable]Could anyone help me understand what is going on? Why the compiler is complaining about those arguments in those functions?
Thank you!!!!
1 Answer
I saved your code as g.c, then issued this command line
gcc -Wall `pkg-config libgvc --cflags --libs` g.c -o EXE -lgvcthat yields
g.c: In function ‘main’:g.c:14:5: error: too few arguments to function ‘agopen’/usr/local/include/graphviz/cgraph.h:266:18: note: declared hereg.c:7:14: warning: unused variable ‘a’ [-Wunused-variable]g.c:6:15: warning: variable ‘e’ set but not used [-Wunused-but-set-variable]then I added the miss parameter
g = agopen("g", Agdirected, 0);and the miss library
gcc -Wall `pkg-config libgvc --cflags --libs` g.c -lgvc -lcgraphnow the code compile and link with just 2 warnings:
g.c: In function ‘main’:g.c:7:14: warning: unused variable ‘a’ [-Wunused-variable]g.c:6:15: warning: variable ‘e’ set but not used [-Wunused-but-set-variable]I think it works because I've built graphviz from source, then pkg-config is up-to-date...The program still need some debug, running it i get:
./a.outThere is no layout engine support for "a.out"Use one of: circo dot fdp neato nop nop1 nop2 osage patchwork sfdp twopiError: Layout was not done. Missing layout plugins? The message is because by default the layout engine use the exe name (i.e. a.out, default of gcc compile-and-link) as layout string...
Agdirected? Is it supposed to be a string? Otherwise you need to declare it. Also, all of the error messages are pretty clear, most of them are because you have to many arguments to some functions. – Some programmer dudeOct 25 '13 at 12:27