C:/cmcintos/defOrgs/examples/vesselCrawler/source/Geom_VesselCrawler.cxx

#ifndef _Geom_VesselCrawler_cxx
#define _Geom_VesselCrawler_cxx
#include "Geom_VesselCrawler.h"

namespace mial
{

        template <class dType, int nDims,class MType, class VType >
        void 
                Geom_VesselCrawler<dType, nDims,MType,VType>
                ::updateMatrixConnections()
        {
                //Matrix connections can not be updated from the mesh,
                //as the mesh only has triangles and thus does not hold all connections

                //TODO throw an exception saying this function cannot be called
        }
        template < class dType, int nDims, class MType, class VType >
        typename Geom_VesselCrawler<dType,nDims,MType,VType>::VectorType
                Geom_VesselCrawler< dType, nDims, MType, VType >
                ::getActiveSprings(int cls)
        {
                //Active springs are considered to be those connected to
                //the front-most layer
                int cnt =0;
                for(int a=0;a<numConnections;a++)
                        if( (connectionClass(a)==cls || cls==-1) && (connectionLayers(a)==numLayers ) )
                                cnt++;

                VectorType springs(cnt);
                cnt=0;
                for(int a=0;a<numConnections;a++)
                        if( (connectionClass(a)==cls || cls==-1) && (connectionLayers(a)==numLayers ) )
                                springs(cnt) = a;

                return springs;
        }

        template < class dType, int nDims, class MType, class VType >
        typename Geom_VesselCrawler<dType,nDims,MType,VType>::VectorType
                Geom_VesselCrawler< dType, nDims, MType, VType >
		::getActiveNodes(int cls)
        {
                //Active nodes are defined as those on the frontmost layer
                //i.e the numNodesPerLayer'th last elements
                int cnt =0;
                for(int a=numNodes;a>numNodes-numNodesPerLayer;a--)
                        if(nodeClass(a)==cls || cls==-1 )
                                cnt++;

                VectorType nodes(cnt);
                cnt=0;
                for(int a=numNodes;a>numNodes-numNodesPerLayer;a--)
                        if(nodeClass(a)==cls || cls ==-1)
                                nodes(cnt) = a;

                return nodes;
        }

        template < class dType, int nDims, class MType, class VType >
        bool
                Geom_VesselCrawler< dType, nDims, MType, VType >
		::addConnection(int start,int end,int cls, int layer)
        {
                //std::cout << start << " " << end << " " << cls << " " << layer << std::endl;
                MatrixType tmpCon;
                tmpCon.set_size(numConnections+1,2);
                VectorType tmpCls(numConnections+1);
                VectorType tmpLayer(numConnections+1);

                //Copy the information
                //tmpCon.update(mConnections,0,0);
                //tmpCls.update(connectionClass,0);
                //tmpLayer.update(connectionLayers,0);
                
                for(int i=0;i<numConnections;i++)
                {
                        tmpCon(i,0) = mConnections(i,0);
                        tmpCon(i,1) = mConnections(i,1);
                        tmpCls(i) = connectionClass(i);
                        tmpLayer(i) = connectionLayers(i);
                }
                tmpCon(numConnections,0) = start;
                tmpCon(numConnections,1) = end;
                tmpCls(numConnections) = cls;
                tmpLayer(numConnections) = layer;

                numConnections++;

                //mConnections.set_size(numConnections,2);
                mConnections = tmpCon;//Copy back the old information plus the new information
                //connectionClass.set_size(numConnections);
                connectionClass = tmpCls;
                //connectionLayers.set_size(numConnections);
                connectionLayers = tmpLayer;

                topologyChange = true;
                return false;
        }

        template < class dType, int nDims, class MType, class VType >
        bool
                Geom_VesselCrawler< dType, nDims, MType, VType >
		::addLayer(MType nodes, VType classes )
        {
                this->addNodes(nodes,classes);
                numLayers++;

                //%%Note e always refers to the node#
                //%Link the new layer to the previous layer (if one exists)
                if(numLayers>1)
                {
                        for (int e=0; e<numNodesPerLayer;e++)
                        {
                                this->addConnection(e+((numLayers-2)*numNodesPerLayer),e+((numLayers-1)*numNodesPerLayer,3,numLayers));
                                //springsClass(a(1)+1,:) = [numLayers-1 numLayers];
                        }
                        //Stability Springs
                        for(int e=1; e<numNodesPerLayer;e++)
                        {
                                //    a = size(springsNodes);
                                if(e+numNodesPerLayer-1<=1+numNodesPerLayer)
                                {
                                        this->addConnection(e+((numLayers-2)*numNodesPerLayer), numNodes,3,numLayers);
                                        //springsClass(a(1)+1,:) = [-3 numLayers];
                                }
                                else{
                                        this->addConnection(e+((numLayers-2)*numNodesPerLayer), e+((numLayers-1)*numNodesPerLayer)-1,3,numLayers);
                                        //springsClass(a(1)+1,:) = [-3 numLayers];
                                }
                                if(e+((numLayers-1)*numNodesPerLayer) < numNodes)
                                {
                                        this->addConnection(e+((numLayers-2)*numNodesPerLayer), e+((numLayers-1)*numNodesPerLayer)+1,3,numLayers);
                                        //springsClass(end+1,:) = [-3 numLayers];
                                }
                                else
                                {
                                        this->addConnection(e+((numLayers-2)*numNodesPerLayer), e+((numLayers-2)*numNodesPerLayer)+2,3,numLayers);
                                        //springsClass(end+1,:) = [-3 numLayers];
                                }
                        }
                }
                //%%%%%%%%%%%%%%%%%%%%%%%%%
                //%%Radial Springs (thickness)
                for(int e=1; e<numNodesPerLayer;e++)
                {
                        this->addConnection(1+(numLayers-1)*numNodesPerLayer, e+((numLayers-1)*numNodesPerLayer),1,numLayers);
                        //springsClass(a(1)+1,:) = [-1 numLayers];
                }
                //%%Boundary Springs
                for(int e = 0; e<numNodesPerLayer-2;e++)
                {
                        this->addConnection(e+1+(numLayers-1)*numNodesPerLayer, e+2+(numLayers-1)*numNodesPerLayer,2,numLayers);
                        //springsClass(a(1)+1,:) = [-2 numLayers];
                }
                this->addConnection(numLayers*numNodesPerLayer,2+(numLayers-1)*numNodesPerLayer,2,numLayers);
                //springsClass(a(1)+1,:) = [-2 numLayers];
                return false;
        }
}//end mial
#endif

---