CylMesh.cpp

Go to the documentation of this file.

/**********************************************************************
** This program is part of 'MOOSE', the
** Messaging Object Oriented Simulation Environment.
**           Copyright (C) 2011 Upinder S. Bhalla. and NCBS
** It is made available under the terms of the
** GNU Lesser General Public License version 2.1
** See the file COPYING.LIB for the full notice.
**********************************************************************/

#include "header.h"
#include "SparseMatrix.h"
#include "ElementValueFinfo.h"
#include "../utility/Vec.h"
#include "Boundary.h"
#include "MeshEntry.h"
// #include "Stencil.h"
#include "ChemCompt.h"
#include "MeshCompt.h"
#include "CubeMesh.h"
#include "CylBase.h"
#include "NeuroNode.h"
// #include "NeuroStencil.h"
#include "NeuroMesh.h"
#include "CylMesh.h"
#include "EndoMesh.h"
#include "../utility/numutil.h"
const Cinfo* CylMesh::initCinfo()
{
        // Field Definitions
        static ElementValueFinfo< CylMesh, double > x0(
            "x0",
            "x coord of one end",
            &CylMesh::setX0,
            &CylMesh::getX0
        );
        static ElementValueFinfo< CylMesh, double > y0(
            "y0",
            "y coord of one end",
            &CylMesh::setY0,
            &CylMesh::getY0
        );
        static ElementValueFinfo< CylMesh, double > z0(
            "z0",
            "z coord of one end",
            &CylMesh::setZ0,
            &CylMesh::getZ0
        );
        static ElementValueFinfo< CylMesh, double > r0(
            "r0",
            "Radius of one end",
            &CylMesh::setR0,
            &CylMesh::getR0
        );
        static ElementValueFinfo< CylMesh, double > x1(
            "x1",
            "x coord of other end",
            &CylMesh::setX1,
            &CylMesh::getX1
        );
        static ElementValueFinfo< CylMesh, double > y1(
            "y1",
            "y coord of other end",
            &CylMesh::setY1,
            &CylMesh::getY1
        );
        static ElementValueFinfo< CylMesh, double > z1(
            "z1",
            "z coord of other end",
            &CylMesh::setZ1,
            &CylMesh::getZ1
        );
        static ElementValueFinfo< CylMesh, double > r1(
            "r1",
            "Radius of other end",
            &CylMesh::setR1,
            &CylMesh::getR1
        );
        static ElementValueFinfo< CylMesh, vector< double > > coords(
            "coords",
            "All the coords as a single vector: x0 y0 z0  x1 y1 z1  r0 r1 diffLength",
            &CylMesh::setCoords,
            &CylMesh::getCoords
        );

        static ElementValueFinfo< CylMesh, double > diffLength(
            "diffLength",
            "Length constant to use for subdivisions"
            "The system will attempt to subdivide using compartments of"
            "length diffLength on average. If the cylinder has different end"
            "diameters r0 and r1, it will scale to smaller lengths"
            "for the smaller diameter end and vice versa."
            "Once the value is set it will recompute diffLength as "
            "totLength/numEntries",
            &CylMesh::setDiffLength,
            &CylMesh::getDiffLength
        );

        static ReadOnlyValueFinfo< CylMesh, unsigned int > numDiffCompts(
            "numDiffCompts",
            "Number of diffusive compartments in model",
            &CylMesh::innerGetNumEntries
        );

        static ReadOnlyValueFinfo< CylMesh, double > totLength(
            "totLength",
            "Total length of cylinder",
            &CylMesh::getTotLength
        );

        // MsgDest Definitions

        // Field Elements

    static Finfo* cylMeshFinfos[] = {
        &x0,            // Value
        &y0,            // Value
        &z0,            // Value
        &r0,            // Value
        &x1,            // Value
        &y1,            // Value
        &z1,            // Value
        &r1,            // Value
        &diffLength,            // Value
        &coords,        // Value
        &numDiffCompts,     // ReadOnlyValue
        &totLength,     // ReadOnlyValue
    };

    static string doc[] =
    {
        "Name", "CylMesh",
        "Author", "Upi Bhalla",
        "Description", "Chemical compartment with cylindrical geometry. "
                "Defaults to a uniform cylinder of radius 1 micron, "
                "length 100 microns, and voxel length 1 micron so there "
                "are 100 voxels in the cylinder. "
                "The cylinder can be given a linear taper, by assigning "
                "different radii r0 and r1 to the two ends. ",
    };
    static Dinfo< CylMesh > dinfo;
    static Cinfo cylMeshCinfo (
        "CylMesh",
        ChemCompt::initCinfo(),
        cylMeshFinfos,
        sizeof( cylMeshFinfos ) / sizeof ( Finfo* ),
        &dinfo,
        doc,
        sizeof(doc)/sizeof(string)
    );

    return &cylMeshCinfo;
}

// Basic class Definitions

static const Cinfo* cylMeshCinfo = CylMesh::initCinfo();

// Class stuff.
CylMesh::CylMesh()
    :
        numEntries_( 100 ),
        useCaps_( 0 ),
        isToroid_( false ),
        x0_( 0.0 ),
        y0_( 0.0 ),
        z0_( 0.0 ),
        x1_( 100.0e-6 ),
        y1_( 0.0 ),
        z1_( 0.0 ),
        r0_( 1.0e-6 ),
        r1_( 1.0e-6 ),
        diffLength_( 1.0e-6 ),
        surfaceGranularity_( 0.1 ),
        totLen_( 100.0e-6 ),
        rSlope_( 0.0 ),
        lenSlope_( 0.0 )
{
    ;
}

CylMesh::~CylMesh()
{
    ;
}

// Field assignment stuff

void CylMesh::updateCoords( const Eref& e, const vector< double >& concs )
{
    double temp = sqrt(
        ( x1_ - x0_ ) * ( x1_ - x0_ ) +
        ( y1_ - y0_ ) * ( y1_ - y0_ ) +
        ( z1_ - z0_ ) * ( z1_ - z0_ )
    );

    if ( doubleEq( temp, 0.0 ) ) {
        cout << "Error: CylMesh::updateCoords:\n"
        "total length of compartment = 0 with these parameters\n";
        return;
    }
    totLen_ = temp;


    temp = totLen_ / diffLength_;
    if ( temp < 1.0 ) {
        diffLength_ = totLen_;
        numEntries_ = 1;
    } else {
        numEntries_ = static_cast< unsigned int >( round ( temp ) );
        diffLength_ = totLen_ / numEntries_;
    }
    rSlope_ = ( r1_ - r0_ ) / numEntries_;
    lenSlope_ = diffLength_ * rSlope_ * 2 / ( r0_ + r1_ );

    // dx2_[0] = diffLength_;
    // dx2_[1] = diffLength_;
    buildStencil();
    setChildConcs( e, concs, 0 );
}

void CylMesh::setX0( const Eref& e, double v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    x0_ = v;
    updateCoords( e, childConcs );
}

double CylMesh::getX0( const Eref& e ) const
{
    return x0_;
}

void CylMesh::setY0( const Eref& e, double v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    y0_ = v;
    updateCoords( e, childConcs );
}

double CylMesh::getY0( const Eref& e ) const
{
    return y0_;
}

void CylMesh::setZ0( const Eref& e, double v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    z0_ = v;
    updateCoords( e, childConcs );
}

double CylMesh::getZ0( const Eref& e ) const
{
    return z0_;
}

void CylMesh::setR0( const Eref& e, double v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    r0_ = v;
    updateCoords( e, childConcs );
}

double CylMesh::getR0( const Eref& e ) const
{
    return r0_;
}


void CylMesh::setX1( const Eref& e, double v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    x1_ = v;
    updateCoords( e, childConcs );
}

double CylMesh::getX1( const Eref& e ) const
{
    return x1_;
}

void CylMesh::setY1( const Eref& e, double v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    y1_ = v;
    updateCoords( e, childConcs );
}

double CylMesh::getY1( const Eref& e ) const
{
    return y1_;
}

void CylMesh::setZ1( const Eref& e, double v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    z1_ = v;
    updateCoords( e, childConcs );
}

double CylMesh::getZ1( const Eref& e ) const
{
    return z1_;
}

void CylMesh::setR1( const Eref& e, double v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    r1_ = v;
    updateCoords( e, childConcs );
}

double CylMesh::getR1( const Eref& e ) const
{
    return r1_;
}

void CylMesh::innerSetCoords( const Eref& e, const vector< double >& v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    x0_ = v[0];
    y0_ = v[1];
    z0_ = v[2];

    x1_ = v[3];
    y1_ = v[4];
    z1_ = v[5];

    r0_ = v[6];
    r1_ = v[7];

    diffLength_ = v[8];

    updateCoords( e, childConcs );
}

void CylMesh::setCoords( const Eref& e, vector< double > v )
{
    if ( v.size() < 9 ) {
        cout << "CylMesh::setCoords: Warning: size of argument vec should be >= 9, was " << v.size() << endl;
    }
    innerSetCoords( e, v );
    transmitChange( e );
}

vector< double > CylMesh::getCoords( const Eref& e ) const
{
    vector< double > ret( 9 );

    ret[0] = x0_;
    ret[1] = y0_;
    ret[2] = z0_;

    ret[3] = x1_;
    ret[4] = y1_;
    ret[5] = z1_;

    ret[6] = r0_;
    ret[7] = r1_;

    ret[8] = diffLength_;
    return ret;
}


void CylMesh::setDiffLength( const Eref& e, double v )
{
    vector< double > childConcs;
    getChildConcs( e, childConcs );
    diffLength_ = v;
    updateCoords( e, childConcs );
}

double CylMesh::getDiffLength( const Eref& e ) const
{
    return diffLength_;
}


double CylMesh::getTotLength() const
{
    return totLen_;
}

unsigned int CylMesh::innerGetDimensions() const
{
    return 3;
}

// FieldElement assignment stuff for MeshEntries

unsigned int CylMesh::getMeshType( unsigned int fid ) const
{
    if ( !isToroid_ && useCaps_ && ( fid == 0 || fid == numEntries_ - 1 ) )
        return SPHERE_SHELL_SEG;

    return CYL;
}

unsigned int CylMesh::getMeshDimensions( unsigned int fid ) const
{
    return 3;
}

double CylMesh::getMeshEntryVolume( unsigned int fid ) const
{
    double len0 = diffLength_ * 2 * r0_ / ( r0_ + r1_ );

    double ri = r0_ + (fid + 0.5) * rSlope_;
    double leni = len0 + ( fid + 0.5 ) * lenSlope_;

    return leni * ri * ri * PI;
}

vector< double > CylMesh::getCoordinates( unsigned int fid ) const
{
    vector< double > ret(10, 0.0);
    double len0 = diffLength_ * 2 * r0_ / ( r0_ + r1_ );
    // double len0 = diffLength_ * 2 * ( r0_ + rSlope_ / 0.5) / ( r0_ + r1_ );
    double lenStart = len0 + lenSlope_ * 0.5;

    double axialStart = fid * lenStart + ( ( fid * (fid - 1 ) )/2 ) * lenSlope_;
        // fid * totLen_/numEntries_ + (fid - frac ) * lenSlope_;
    double axialEnd =
        (fid + 1) * lenStart + ( ( fid * (fid + 1 ) )/2 ) * lenSlope_;
        // (fid + 1) * totLen_/numEntries_ + (fid - frac + 1.0) * lenSlope_;

    ret[0] = x0_ + (x1_ - x0_ ) * axialStart/totLen_;
    ret[1] = y0_ + (y1_ - y0_ ) * axialStart/totLen_;
    ret[2] = z0_ + (z1_ - z0_ ) * axialStart/totLen_;

    ret[3] = x0_ + (x1_ - x0_ ) * axialEnd/totLen_;
    ret[4] = y0_ + (y1_ - y0_ ) * axialEnd/totLen_;
    ret[5] = z0_ + (z1_ - z0_ ) * axialEnd/totLen_;

    ret[6] = r0_ + fid * rSlope_;
    ret[7] = r0_ + (fid + 1.0) * rSlope_;

    ret[8] = 0;
    ret[9] = 0;

    return ret;
}

vector< double > CylMesh::getDiffusionArea( unsigned int fid ) const
{
    if ( numEntries_ <= 1 )
        return vector< double >( 0 );

    double rlow = r0_ + fid * rSlope_;
    double rhigh = r0_ + (fid + 1.0) * rSlope_;

    if ( fid == 0 ) {
        if ( isToroid_ ) {
            vector < double > ret( 2 );
            ret[0] = rlow * rlow * PI;
            ret[1] = rhigh * rhigh * PI;
            return ret;
        } else {
            return vector < double >( 1, rhigh * rhigh * PI );
        }
    }

    if ( fid == (numEntries_ - 1 ) ) {
        if ( isToroid_ ) {
            vector < double > ret( 2 );
            ret[0] = rlow * rlow * PI;
            ret[1] = r0_ * r0_ * PI; // Wrapping around
            return ret;
        } else {
            return vector < double >( 1, rlow * rlow * PI );
        }
    }
    vector< double > ret( 2 );
    ret[0] = rlow * rlow * PI;
    ret[1] = rhigh * rhigh * PI;
    return ret;
}

vector< double > CylMesh::getDiffusionScaling( unsigned int fid ) const
{
    if ( numEntries_ <= 1 )
        return vector< double >( 0 );

    if ( !isToroid_ && ( fid == 0 || fid == (numEntries_ - 1) ) )
        return vector< double >( 1, 1.0 );

    return vector< double >( 2, 1.0 );
}

double CylMesh::extendedMeshEntryVolume( unsigned int fid ) const
{
    if ( fid < numEntries_ ) {
        return getMeshEntryVolume( fid );
    } else {
        return MeshCompt::extendedMeshEntryVolume( fid - numEntries_ );
    }
}

// Dest funcsl

void CylMesh::innerHandleRequestMeshStats( const Eref& e,
        const SrcFinfo2< unsigned int, vector< double > >* meshStatsFinfo
    )
{
    vector< double > ret( vGetEntireVolume() / numEntries_ ,1 );
    meshStatsFinfo->send( e, 1, ret );
}

void CylMesh::innerHandleNodeInfo(
            const Eref& e,
            unsigned int numNodes, unsigned int numThreads )
{
        /*
    unsigned int numEntries = numEntries_;
    vector< double > vols( numEntries, volume_ / numEntries );
    vector< unsigned int > localEntries( numEntries );
    vector< vector< unsigned int > > outgoingEntries;
    vector< vector< unsigned int > > incomingEntries;
    */
    /*
    double oldvol = getMeshEntryVolume( 0 );
    meshSplit()->send( e,
        oldvol,
        vols, localEntries,
        outgoingEntries, incomingEntries );
        */
}

unsigned int CylMesh::innerGetNumEntries() const
{
    return numEntries_;
}

void CylMesh::innerSetNumEntries( unsigned int n )
{
    static const unsigned int WayTooLarge = 1000000;
    if ( n == 0 || n > WayTooLarge ) {
        cout << "Warning: CylMesh::innerSetNumEntries( " << n <<
        " ): out of range\n";
        return;
    }
    assert( n > 0 );
    numEntries_ = n;
    diffLength_ = totLen_ / n;
    rSlope_ = ( r1_ - r0_ ) / numEntries_;
    lenSlope_ = diffLength_ * rSlope_ * 2 / ( r0_ + r1_ );

    buildStencil();
}


void CylMesh::innerBuildDefaultMesh( const Eref& e,
    double volume, unsigned int numEntries )
{
    double r = pow( ( volume / ( PI * 2 ) ), 1.0 / 3 );
    vector< double > coords( 9, 0 );
    coords[3] = 2 * r;
    coords[6] = r;
    coords[7] = r;
    coords[8] = 2 * r / numEntries;
    setCoords( e, coords );
}

vector< unsigned int > CylMesh::getParentVoxel() const
{
    vector< unsigned int > ret( numEntries_ );
    if ( numEntries_ > 0 )
        ret[0] = static_cast< unsigned int >( -1 );
    for (unsigned int i = 1; i < numEntries_; ++i )
        ret[i] = i-1;

    return ret;
}

const vector< double >& CylMesh::vGetVoxelVolume() const
{
    static vector< double > vol;
    vol.resize( numEntries_ );
    for ( unsigned int i = 0; i < numEntries_; ++i )
        vol[i] = getMeshEntryVolume( i );
    return vol;
}

const vector< double >& CylMesh::vGetVoxelMidpoint() const
{
    static vector< double > midpoint( numEntries_ * 3, 0.0 );
    midpoint.resize( numEntries_ * 3 );
    double dx = ( x1_ - x0_ ) / numEntries_;
    double dy = ( y1_ - y0_ ) / numEntries_;
    double dz = ( z1_ - z0_ ) / numEntries_;
    vector< double >::iterator j = midpoint.begin();
    for ( unsigned int i = 0; i < numEntries_; ++i )
        *j++ = x0_ + dx * i;
    for ( unsigned int i = 0; i < numEntries_; ++i )
        *j++ = y0_ + dy * i;
    for ( unsigned int i = 0; i < numEntries_; ++i )
        *j++ = z0_ + dz * i;

    return midpoint;
}

const vector< double >& CylMesh::getVoxelArea() const
{
    static vector< double > area;
    area.resize( numEntries_ );
    for ( unsigned int i = 0; i < numEntries_; ++i ) {
        double frac = ( 0.5 + static_cast< double >( i ) ) /
            static_cast< double >( numEntries_ );
        double r = r0_ * ( 1.0 - frac ) + r1_ * frac;
        area[i] = r * r * PI;
    }
    return area;
}

const vector< double >& CylMesh::getVoxelLength() const
{
    static vector< double > length;
    length.assign( numEntries_, totLen_ / numEntries_ );
    return length;
}

double CylMesh::vGetEntireVolume() const
{
    double vol = 0.0;
    for ( unsigned int i = 0; i < numEntries_; ++i )
        vol += getMeshEntryVolume( i );
    return vol;
}

bool CylMesh::vSetVolumeNotRates( double volume )
{
    double oldVol = vGetEntireVolume();
    double linScale = pow( volume/oldVol, 1.0 / 3.0 );
    x1_ *= linScale;
    y1_ *= linScale;
    z1_ *= linScale;
    r0_ *= linScale;
    r1_ *= linScale;
    totLen_ *= linScale;
    // Have to scale this so numEntries remains the same.
    diffLength_ = totLen_ / numEntries_;
    return true;
}

// Utility function to transmit any changes to target nodes.

void CylMesh::transmitChange( const Eref& e )
{
        /*
    Id meshEntry( e.id().value() + 1 );
    assert(
        meshEntry.eref().data() == reinterpret_cast< char* >( lookupEntry( 0 ) )
    );
    double oldvol = getMeshEntryVolume( 0 );
    unsigned int totalNumEntries = numEntries_;
    unsigned int localNumEntries = totalNumEntries;
    unsigned int startEntry = 0;
    vector< unsigned int > localIndices( localNumEntries ); // empty
    for ( unsigned int i = 0; i < localNumEntries; ++i )
        localIndices[i] = i;
    vector< double > vols( localNumEntries, volume_ / numEntries_ );
    vector< vector< unsigned int > > outgoingEntries; // [node#][Entry#]
    vector< vector< unsigned int > > incomingEntries; // [node#][Entry#]

    // This function updates the size of the FieldDataHandler for the
    // MeshEntries.
    DataHandler* dh = meshEntry.element()->dataHandler();
    FieldDataHandlerBase* fdh = dynamic_cast< FieldDataHandlerBase* >( dh );
    assert( fdh );
    if ( totalNumEntries > fdh->getMaxFieldEntries() ) {
        fdh->setMaxFieldEntries( localNumEntries );
    }

    // This message tells the Stoich about the new mesh, and also about
    // how it communicates with other nodes.
    meshSplit()->fastSend( e,
        oldvol,
        vols, localIndices,
        outgoingEntries, incomingEntries );

    // This func goes down to the MeshEntry to tell all the pools and
    // Reacs to deal with the new mesh. They then update the stoich.
    lookupEntry( 0 )->triggerRemesh( meshEntry.eref(),
        oldvol, startEntry, localIndices, vols );
        */
}

// Utility function to set up Stencil for diffusion
void CylMesh::buildStencil()
{
    setStencilSize( numEntries_, numEntries_ );
    for ( unsigned int i = 0; i < numEntries_; ++i ) {
        double rLow = r0_ + i * rSlope_;
        double rHigh = r0_ + (i + 1.0) * rSlope_;
        double aLow = rLow * rLow * PI;
        double aHigh = rHigh * rHigh * PI;
        vector< double > entry;
        vector< unsigned int > colIndex;
        if ( i == 0 ) {
            colIndex.push_back( 1 );
            entry.push_back( aHigh / diffLength_ );
            if ( isToroid_ ) {
                colIndex.push_back( numEntries_ - 1 );
                entry.push_back( aLow / diffLength_ );
            }
        } else if ( i == numEntries_ - 1 ) {
            if ( isToroid_ ) {
                colIndex.push_back( 0 );
                if ( r0_ < r1_ )
                    entry.push_back( r0_ * r0_ * PI / diffLength_ );
                else
                    entry.push_back( r1_ * r1_ * PI / diffLength_ );
            }
            colIndex.push_back( numEntries_ - 2 );
            entry.push_back( aLow / diffLength_ );
        } else { // Mostly it is in the middle.
            colIndex.push_back( i - 1 );
            entry.push_back( aLow / diffLength_ );
            colIndex.push_back( i + 1 );
            entry.push_back( aHigh / diffLength_ );
        }
        addRow( i, entry, colIndex );
    }
    innerResetStencil();
}

// Utility function for junctions

void CylMesh::matchMeshEntries( const ChemCompt* other,
       vector< VoxelJunction >& ret ) const
{
    // This is seriously ugly, and what virtual funcs were meant to handle.
    // Dirty hack for now.
    const CylMesh* cyl = dynamic_cast< const CylMesh* >( other );
    if ( cyl ) {
        matchCylMeshEntries( cyl, ret );
        return;
    }
    const EndoMesh* em = dynamic_cast< const EndoMesh* >( other );
    if ( em )
    {
        em->matchMeshEntries( this, ret );
        flipRet( ret );
        return;
    }
    const CubeMesh* cube = dynamic_cast< const CubeMesh* >( other );
    if ( cube ) {
        matchCubeMeshEntries( cube, ret );
        return;
    }
    const NeuroMesh* nm = dynamic_cast< const NeuroMesh* >( other );
    if ( nm ) {
        matchNeuroMeshEntries( nm, ret );
        return;
    }
    cout << "Warning:CylMesh::matchMeshEntries: "  << 
            " unknown mesh type\n";
}

// Look for end-to-end diffusion, not sideways for now.
// Very straightforward but tedious because of the different permutations.
// Could readily add cylinder side to other end.
void CylMesh::matchCylMeshEntries( const CylMesh* other,
vector< VoxelJunction >& ret ) const
{
    const double EPSILON = 1e-3;
    ret.clear();
    // Should really estimate the distance from the centre of the smaller
    // cylinder cap disk to the plane of the larger disk, provided it is
    // within the radius of the disk. The subsequent calculations are the
    // same though.
    double dr1 =    // startOfSelf-to-startOfOther
            distance(x0_ - other->x0_, y0_ - other->y0_, z0_ - other->z0_ );
    double dr2 =  // endOfSelf-to-endOfOther
            distance(x1_ - other->x1_, y1_ - other->y1_, z1_ - other->z1_ );
    double dr3 = // endOfSelf-to-startOfOther
            distance(x1_ - other->x0_, y1_ - other->y0_, z1_ - other->z0_ );
    double dr4 =  // startOfSelf-to-endOfOther
            distance(x0_ - other->x1_, y0_ - other->y1_, z0_ - other->z1_ );

    if ( dr1 <= dr2 && dr1 <= dr3 && dr1 <= dr4 ) {
        if ( ( dr1/totLen_ < EPSILON && dr1/other->totLen_ < EPSILON ) ) {
            double xda;
            if ( r0_ < other->r0_ )
                xda = 2 * r0_ * r0_ * PI / ( diffLength_ + other->diffLength_ );
            else
                xda = 2 * other->r0_ * other->r0_ * PI /
                        ( diffLength_ + other->diffLength_ );
            ret.push_back( VoxelJunction( 0, 0, xda ) );
            ret.back().first = 0;
            ret.back().second = 0;
            ret.back().firstVol = getMeshEntryVolume( 0 );
            ret.back().secondVol = other->getMeshEntryVolume( 0 );
        }
    } else if ( dr2 <= dr3 && dr2 <= dr4 ) {
        if ( ( dr2/totLen_ < EPSILON && dr2/other->totLen_ < EPSILON ) ) {
            double xda;
            if ( r1_ < other->r1_ )
                xda = 2 * r1_ * r1_ * PI / ( diffLength_ + other->diffLength_ );
            else
                xda = 2 * other->r1_ * other->r1_ * PI /
                        ( diffLength_ + other->diffLength_ );
            ret.push_back( VoxelJunction( numEntries_ - 1,
                        other->numEntries_ - 1, xda ) );
            ret.back().first = numEntries_;
            ret.back().second = other->numEntries_ - 1;
            ret.back().firstVol = getMeshEntryVolume( numEntries_ - 1 );
            ret.back().secondVol = other->getMeshEntryVolume( other->numEntries_ - 1 );
        }
    } else if ( dr3 <= dr4 ) {
        if ( ( dr3/totLen_ < EPSILON && dr3/other->totLen_ < EPSILON ) ) {
            double xda;
            if ( r1_ < other->r0_ )
                xda = 2 * r1_ * r1_ * PI / ( diffLength_ + other->diffLength_ );
            else
                xda = 2 * other->r0_ * other->r0_ * PI /
                        ( diffLength_ + other->diffLength_ );
            ret.push_back( VoxelJunction( numEntries_ - 1, 0, xda ) );
            ret.back().first = numEntries_ - 1;
            ret.back().second = 0;
            ret.back().firstVol = getMeshEntryVolume( numEntries_ - 1 );
            ret.back().secondVol = other->getMeshEntryVolume( 0 );
        }
    } else {
        if ( ( dr4/totLen_ < EPSILON && dr4/other->totLen_ < EPSILON ) ) {
            double xda;
            if ( r0_ < other->r1_ )
                xda = 2 * r0_ * r0_ * PI / ( diffLength_ + other->diffLength_ );
            else
                xda = 2 * other->r1_ * other->r1_ * PI /
                        ( diffLength_ + other->diffLength_ );
            ret.push_back( VoxelJunction( 0, other->numEntries_ - 1, xda ));
            ret.back().first = 0;
            ret.back().second = other->numEntries_ - 1;
            ret.back().firstVol = getMeshEntryVolume( 0 );
            ret.back().secondVol = other->getMeshEntryVolume( other->numEntries_ - 1 );
        }
    }
}

// Select grid volume. Ideally the meshes should be comparable.
double CylMesh::selectGridVolume( double h ) const
{
    if ( h > diffLength_ )
        h = diffLength_;
    if ( h > r0_ )
        h = r0_;
    if ( h > r1_ )
        h = r1_;
    h *= surfaceGranularity_;
    unsigned int num = ceil( diffLength_ / h );
    h = diffLength_ / num;

    return h;
}

void fillPointsOnCircle(
                const Vec& u, const Vec& v, const Vec& q,
                double h, double r, vector< double >& area,
                const CubeMesh* other
                )
{
    // fine-tune the h spacing so it is integral around circle.
    // This will cause small errors in area estimate but they will
    // be anisotropic. The alternative will have large errors toward
    // 360 degrees, but not elsewhere.
    unsigned int numAngle = floor( 2.0 * PI * r / h + 0.5 );
    assert( numAngle > 0 );
    double dtheta = 2.0 * PI / numAngle;
    double dArea = h * dtheta * r;
    // March along points on surface of circle centred at q.
    for ( unsigned int j = 0; j < numAngle; ++j ) {
        double theta = j * dtheta;
        double c = cos( theta );
        double s = sin( theta );
        double p0 = q.a0() + r * ( u.a0() * c + v.a0() * s );
        double p1 = q.a1() + r * ( u.a1() * c + v.a1() * s );
        double p2 = q.a2() + r * ( u.a2() * c + v.a2() * s );
        unsigned int index = other->spaceToIndex( p0, p1, p2 );
        if ( index != CubeMesh::EMPTY )
            area[index] += dArea;
    }
}

void CylMesh::matchCubeMeshEntries( const CubeMesh* other,
vector< VoxelJunction >& ret ) const
{
    const double EPSILON = 1e-18;
    Vec a( x1_ - x0_, y1_ - y0_, z1_ - z0_ );
    Vec u;
    Vec v;
    a.orthogonalAxes( u, v );

    double h = selectGridVolume( other->getDx() );

    unsigned int num = floor( 0.1 + diffLength_ / h );
    // March along axis of cylinder.
    // q is the location of the point along axis.
    for ( unsigned int i = 0; i < numEntries_; ++i ) {
        vector< double >area( other->getNumEntries(), 0.0 );
        for ( unsigned int j = 0; j < num; ++j ) {
            unsigned int m = i * num + j;
            double frac = ( m * h + h/2.0 ) / totLen_;
            double q0 = x0_ + a.a0() * frac;
            double q1 = y0_ + a.a1() * frac;
            double q2 = z0_ + a.a2() * frac;
            // get radius of cylinder at this point.
            double r = r0_ + ( m * h + h / 2.0 ) * rSlope_;
            fillPointsOnCircle( u, v, Vec( q0, q1, q2 ),
                        h, r, area, other );
            }
        // Go through all cubeMesh entries and compute diffusion
        // cross-section. Assume this is through a membrane, so the
        // only factor relevant is area. Not the distance.
        for ( unsigned int k = 0; k < area.size(); ++k ) {
            if ( area[k] > EPSILON ) {
                ret.push_back( VoxelJunction( i, k, area[k] ) );
            }
        }
    }
}

void CylMesh::matchNeuroMeshEntries( const NeuroMesh* other,
vector< VoxelJunction >& ret ) const
{
}

void CylMesh::indexToSpace( unsigned int index,
            double& x, double& y, double& z ) const
{
    if ( index < numEntries_ ) {
        double k = ( index + 0.5 ) / static_cast< double >( numEntries_ );
        x = ( x1_ - x0_ ) * k + x0_;
        y = ( y1_ - y0_ ) * k + y0_;
        z = ( z1_ - z0_ ) * k + z0_;
    }
}

static double dotprd ( double x0, double y0, double z0,
                double x1, double y1, double z1 )
{
        return x0 * x1 + y0 * y1 + z0 * z1;
}


// this is the function that does the actual calculation.
double CylMesh::nearest( double x, double y, double z,
                double& linePos, double& r ) const
{
    // Consider r0 = x0,y0,z0 and r1 = x1, y1, z1, and r = x,y,z.
    // Fraction along cylinder = k
    //
    // Then, point p along line from r0 to r1 is
    // p = k( r0-r1) + r1.
    //
    // Solving,
    // k = (r0 - r1).(r - r1) / (|r0-r1|^2)
    //

    double dist = distance( x1_ - x0_, y1_ - y0_, z1_ - z0_ );
    double k = dotprd(
        x1_ - x0_, y1_ - y0_, z1_ - z0_,
        x - x0_, y - y0_, z - z0_ ) / ( dist * dist );

    // x2, y2, z2 are the coords of the nearest point.
    double x2 = k * (x1_ - x0_) + x0_;
    double y2 = k * (y1_ - y0_) + y0_;
    double z2 = k * (z1_ - z0_) + z0_;

    double ret = distance( x - x2, y - y2, z - z2 );
    linePos = k;
    r = r0_ + k * numEntries_ * rSlope_;
    return ret;
}

// This function returns the index.
double CylMesh::nearest( double x, double y, double z,
                unsigned int& index ) const
{
    double k = 0.0;
    double r;
    double ret = nearest( x, y, z, k, r );
    if ( k < 0.0 ) {
        ret = -ret;
        index = 0;
    } else if ( k > 1.0 ) {
        ret = -ret;
        index = numEntries_ - 1;
    } else { // Inside length of cylinder, now is it inside radius?
        index = k * numEntries_;
        double ri = r0_ + (index + 0.5) * rSlope_;
        if ( ret > ri )
            ret = -ret;
    }
    return ret;
}


/*
bool isOnSurface( double x, double y, double z,
                    double dx, double dy, double dz,
                    unsigned int &index, double& adx )
{
    double len = distance( x1_ - x0_, y1_ - y0_, z1_ - z0_ );
    double k = dotprd(
        x1_ - x0_, y1_ - y0_, z1_ - z0_,
        x - x0_, y - y0_, z - z0_ ) / len;

    // x2, y2, z2 are the coords of the nearest point.
    double x2 = k * (x1_ - x0_) + x0_;
    double y2 = k * (y1_ - y0_) + y0_;
    double z2 = k * (z1_ - z0_) + z0_;

    double ret = distance( x - x2, y - y2, z - z2 );

    double cubeRange = sqrt(dx*dx + dy*dy + dz*dz);

    // Now we check if the distance is definitely too far off for the
    // passed in point
    if ( k < -dx/2 || k > len + dx/2 ) // past the end.
        return false;

    double ri = k * rSlope_; // local cylinder radius.

    if ( ret > ri + cubeRange || ret < ri - cubeRange )
        return false;

    // OK, now we need to find the plane of intersection of the cylinder
    // with the cuboid. To make it easier, assume it is flat. We already
    // know the vector from the middle of the cuboid to the nearest
    // cylinder point. Treat it as the normal to the intersection plane.
    // We need: : is the plane inside the cube?
    // What is the area of the plane till its intersection with the cube?

}

*/