testMesh.cpp File Reference

#include "header.h"
#include "SparseMatrix.h"
#include "../shell/Shell.h"
#include "Boundary.h"
#include "MeshEntry.h"
#include "ChemCompt.h"
#include "MeshCompt.h"
#include "CubeMesh.h"
#include "CylBase.h"
#include "NeuroNode.h"
#include "NeuroMesh.h"
#include "../utility/Vec.h"
#include "CylMesh.h"
#include "SpineEntry.h"
#include "SpineMesh.h"
#include "PsdMesh.h"

Include dependency graph for testMesh.cpp:

Go to the source code of this file.

Typedefs
typedef pair< unsigned int, unsigned int >	PII

Functions
void	checkAbut (const vector< PII > &intersect, unsigned int ix, unsigned int iy, unsigned int iz, unsigned int nx, unsigned int ny, unsigned int nz, unsigned int meshIndex, vector< VoxelJunction > &ret)
Id	makeCompt (Id parentCompt, Id parentObj, string name, double len, double dia, double theta)
void	setIntersectVoxel (vector< PII > &intersect, unsigned int ix, unsigned int iy, unsigned int iz, unsigned int nx, unsigned int ny, unsigned int nz, unsigned int meshIndex)
void	testCubeMesh ()
void	testCubeMeshExtendStencil ()
void	testCubeMeshFillThreeDimSurface ()
void	testCubeMeshFillTwoDimSurface ()
void	testCubeMeshJunctionDiffSizeMesh ()
void	testCubeMeshJunctionThreeDimSurface ()
void	testCubeMeshJunctionTwoDimSurface ()
void	testCubeMeshMultiJunctionTwoD ()
void	testIntersectVoxel ()
void	testMesh ()
void	testReMesh ()
void	testVec ()
void	testVolScaling ()

Variables
static const unsigned int	ABUT = ~2
static const unsigned int	EMPTY = ~0
static const unsigned int	MULTI = ~3
static const unsigned int	SURFACE = ~1

Typedef Documentation

typedef pair< unsigned int, unsigned int > PII

Definition at line 1304 of file testMesh.cpp.

Function Documentation

void checkAbut	(	const vector< PII > &	intersect,
		unsigned int	ix,
		unsigned int	iy,
		unsigned int	iz,
		unsigned int	nx,
		unsigned int	ny,
		unsigned int	nz,
		unsigned int	meshIndex,
		vector< VoxelJunction > &	ret
	)

checkAbut checks the intersect vector for the current position ix, iy, iz, to determine how many diffusion terms to extract. It then puts each of the extracted terms into the ret vector. There is a minor efficiency for one and two diffusion terms as they are encoded within the intersect vector. Higher-order surface alignments require an in-line scan of neighboring voxels. In all casesl the function inserts a flag indicating surface direction into the diffScale field of the VoxelJunction. 0 = x; 1 = y; 2 = z.

Definition at line 1348 of file CubeMesh.cpp.

References CubeMesh::ABUTX, CubeMesh::ABUTY, CubeMesh::ABUTZ, CubeMesh::EMPTY, CubeMesh::MULTI, and CubeMesh::SURFACE.

Referenced by CubeMesh::matchCubeMeshEntries(), and testIntersectVoxel().

{
    unsigned int index = ( iz * ny + iy ) * nx + ix;
    unsigned int localFlag = intersect[index].second;

    if ( localFlag == CubeMesh::EMPTY || localFlag == CubeMesh::SURFACE )
            return; // Nothing to put into the ret vector
    if ( localFlag == CubeMesh::ABUTX ) {
        ret.push_back(
                VoxelJunction( intersect[index].first, meshIndex, 0 ));
    } else if ( localFlag == CubeMesh::ABUTY ) {
        ret.push_back(
                VoxelJunction( intersect[index].first, meshIndex, 1 ));
    } else if ( localFlag == CubeMesh::ABUTZ ) {
        ret.push_back(
                VoxelJunction( intersect[index].first, meshIndex, 2 ));
    } else if ( localFlag == CubeMesh::MULTI ) { // go through all 6 cases.
        if ( ix > 0 ) {
            index = ( iz * ny + iy ) * nx + ix - 1;
            if ( intersect[index].second == CubeMesh::SURFACE )
                ret.push_back(
                    VoxelJunction( intersect[index].first, meshIndex, 0 ));
        }
        if ( ix + 1 < nx ) {
            index = ( iz * ny + iy ) * nx + ix + 1;
            if ( intersect[index].second == CubeMesh::SURFACE )
                ret.push_back(
                    VoxelJunction( intersect[index].first, meshIndex, 0 ));
        }
        if ( iy > 0 ) {
            index = ( iz * ny + iy -1 ) * nx + ix;
            if ( intersect[index].second == CubeMesh::SURFACE )
                ret.push_back(
                    VoxelJunction( intersect[index].first, meshIndex, 1 ));
        }
        if ( iy + 1 < ny ) {
            index = ( iz * ny + iy + 1 ) * nx + ix;
            if ( intersect[index].second == CubeMesh::SURFACE )
                ret.push_back(
                    VoxelJunction( intersect[index].first, meshIndex, 1 ));
        }
        if ( iz > 0 ) {
            index = ( (iz-1) * ny + iy ) * nx + ix;
            if ( intersect[index].second == CubeMesh::SURFACE )
                ret.push_back(
                    VoxelJunction( intersect[index].first, meshIndex, 2 ));
        }
        if ( iz + 1 < nz ) {
            index = ( (iz+1) * ny + iy ) * nx + ix;
            if ( intersect[index].second == CubeMesh::SURFACE )
                ret.push_back(
                    VoxelJunction( intersect[index].first, meshIndex, 2 ));
        }
    }
}

Here is the caller graph for this function:

Id makeCompt	(	Id	parentCompt,
		Id	parentObj,
		string	name,
		double	len,
		double	dia,
		double	theta
	)

theta in degrees len and dia in metres as usual

Definition at line 855 of file testMesh.cpp.

References Eref::data(), Shell::doAddMsg(), Shell::doCreate(), Id::eref(), Field< A >::get(), PI, and Field< A >::set().

Referenced by testCellDiffn(), and testSmallCellDiffn().

{
    Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
    Id ret = shell->doCreate( "Compartment", parentObj, name, 1 );
    double pax = 0;
    double pay = 0;
    if ( parentCompt != Id() ) {
        pax = Field< double >::get( parentCompt, "x" );
        pay = Field< double >::get( parentCompt, "y" );
        shell->doAddMsg( "Single", parentCompt, "raxial", ret, "axial" );
    }
    Field< double >::set( ret, "x0", pax );
    Field< double >::set( ret, "y0", pay );
    Field< double >::set( ret, "z0", 0.0 );
    double x = pax + len * cos( theta * PI / 180.0 );
    double y = pay + len * sin( theta * PI / 180.0 );
    Field< double >::set( ret, "x", x );
    Field< double >::set( ret, "y", y );
    Field< double >::set( ret, "z", 0.0 );
    Field< double >::set( ret, "diameter", dia );
    Field< double >::set( ret, "length", len );

    return ret;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void setIntersectVoxel	(	vector< PII > &	intersect,
		unsigned int	ix,
		unsigned int	iy,
		unsigned int	iz,
		unsigned int	nx,
		unsigned int	ny,
		unsigned int	nz,
		unsigned int	meshIndex
	)

Definition at line 1304 of file CubeMesh.cpp.

References CubeMesh::ABUTX, CubeMesh::ABUTY, CubeMesh::ABUTZ, setAbut(), and CubeMesh::SURFACE.

Referenced by CubeMesh::assignVoxels(), and testIntersectVoxel().

{
    assert( ix < nx && iy < ny && iz < nz );
    unsigned int index = ( iz * ny + iy ) * nx + ix;
    assert( index < intersect.size() );
    intersect[index] = PII( meshIndex, CubeMesh::SURFACE );
    if ( ix > 0 )
        setAbut( intersect[ (iz*ny + iy) * nx + ix-1 ], meshIndex,
                        CubeMesh::ABUTX );
    if ( ix + 1 < nx )
        setAbut( intersect[ (iz*ny + iy) * nx + ix+1 ], meshIndex,
                        CubeMesh::ABUTX );

    if ( iy > 0 )
        setAbut( intersect[ ( iz*ny + iy-1 ) * nx + ix ], meshIndex,
                        CubeMesh::ABUTY);
    if ( iy + 1 < ny )
        setAbut( intersect[ ( iz*ny + iy+1 ) * nx + ix ], meshIndex,
                        CubeMesh::ABUTY);

    if ( iz > 0 )
        setAbut( intersect[ ( (iz-1)*ny + iy ) * nx + ix ], meshIndex,
                        CubeMesh::ABUTZ);
    if ( iz + 1 < nz )
        setAbut( intersect[ ( (iz+1)*ny + iy ) * nx + ix ], meshIndex,
                        CubeMesh::ABUTZ);
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testCubeMesh

(

)

Low-level tests for the CubeMesh object: No MOOSE calls involved.

Definition at line 600 of file testMesh.cpp.

References CUBOID, doubleEq(), CubeMesh::getCoords(), ChemCompt::getDimensions(), CubeMesh::getDx(), CubeMesh::getDy(), CubeMesh::getDz(), CubeMesh::getMeshDimensions(), CubeMesh::getMeshType(), MeshCompt::getNeighbors(), CubeMesh::getNx(), CubeMesh::getNy(), CubeMesh::getNz(), CubeMesh::getX0(), CubeMesh::getX1(), CubeMesh::getY0(), CubeMesh::getY1(), CubeMesh::getZ0(), CubeMesh::getZ1(), CubeMesh::innerGetNumEntries(), CubeMesh::innerSetCoords(), CubeMesh::setPreserveNumEntries(), CubeMesh::setX0(), CubeMesh::setX1(), CubeMesh::setY0(), CubeMesh::setY1(), CubeMesh::setZ0(), and CubeMesh::setZ1().

Referenced by testMesh().

{
    CubeMesh cm;
    cm.setPreserveNumEntries( 0 );
    assert( cm.getMeshType( 0 ) == CUBOID );
    assert( cm.getMeshDimensions( 0 ) == 3 );
    assert( cm.getDimensions() == 3 );

    vector< double > coords( 9 );
    coords[0] = 0; // X0
    coords[1] = 0; // Y0
    coords[2] = 0; // Z0

    coords[3] = 2; // X1
    coords[4] = 4; // Y1
    coords[5] = 8; // Z1

    coords[6] = 1; // DX
    coords[7] = 1; // DY
    coords[8] = 1; // DZ

    cm.innerSetCoords( coords );

    vector< unsigned int > neighbors = cm.getNeighbors( 0 );
    assert( neighbors.size() == 3 );
    assert( neighbors[0] = 1 );
    assert( neighbors[0] = 2 );
    assert( neighbors[0] = 8 );

    assert( cm.innerGetNumEntries() == 64 );
    assert( doubleEq( cm.getX0(), 0 ) );
    assert( doubleEq( cm.getY0(), 0 ) );
    assert( doubleEq( cm.getZ0(), 0 ) );

    assert( doubleEq( cm.getX1(), 2 ) );
    assert( doubleEq( cm.getY1(), 4 ) );
    assert( doubleEq( cm.getZ1(), 8 ) );

    assert( doubleEq( cm.getDx(), 1 ) );
    assert( doubleEq( cm.getDy(), 1 ) );
    assert( doubleEq( cm.getDz(), 1 ) );

    assert( cm.getNx() == 2 );
    assert( cm.getNy() == 4 );
    assert( cm.getNz() == 8 );

    cm.setX0( 1 );
    cm.setY0( 2 );
    cm.setZ0( 4 );

    cm.setX1( 5 );
    cm.setY1( 6 );
    cm.setZ1( 8 );

    vector< double > temp = cm.getCoords( Id().eref() );
    assert( temp.size() == 9 );
    assert( doubleEq( temp[0], 1 ) );
    assert( doubleEq( temp[1], 2 ) );
    assert( doubleEq( temp[2], 4 ) );
    assert( doubleEq( temp[3], 5 ) );
    assert( doubleEq( temp[4], 6 ) );
    assert( doubleEq( temp[5], 8 ) );
    assert( doubleEq( temp[6], 1 ) );
    assert( doubleEq( temp[7], 1 ) );
    assert( doubleEq( temp[8], 1 ) );
    assert( cm.innerGetNumEntries() == 64 );
    assert( cm.getNx() == 4 );
    assert( cm.getNy() == 4 );
    assert( cm.getNz() == 4 );

    neighbors = cm.getNeighbors( 0 );
    assert( neighbors.size() == 3 );
    assert( neighbors[0] == 1 );
    assert( neighbors[1] == 4 );
    assert( neighbors[2] = 16 );

    neighbors = cm.getNeighbors( 63 );
    assert( neighbors.size() == 3 );
    assert( neighbors[0] == 47 );
    assert( neighbors[1] == 59 );
    assert( neighbors[2] == 62 );

    neighbors = cm.getNeighbors( 2 );
    assert( neighbors.size() == 4 );
    assert( neighbors[0] == 1 );
    assert( neighbors[1] == 3 );
    assert( neighbors[2] == 6 );
    assert( neighbors[3] == 18 );

    neighbors = cm.getNeighbors( 6 );
    assert( neighbors.size() == 5 );
    assert( neighbors[0] == 2 );
    assert( neighbors[1] == 5 );
    assert( neighbors[2] == 7 );
    assert( neighbors[3] == 10 );
    assert( neighbors[4] == 22 );

    neighbors = cm.getNeighbors( 22 );
    assert( neighbors.size() == 6 );
    assert( neighbors[0] == 6 );
    assert( neighbors[1] == 18 );
    assert( neighbors[2] == 21 );
    assert( neighbors[3] == 23 );
    assert( neighbors[4] == 26 );
    assert( neighbors[5] == 38 );

    cm.setPreserveNumEntries( 1 );
    assert( cm.getNx() == 4 );
    assert( cm.getNy() == 4 );
    assert( cm.getNz() == 4 );
    assert( doubleEq( cm.getDx(), 1.0 ) );
    assert( doubleEq( cm.getDy(), 1.0 ) );
    assert( doubleEq( cm.getDz(), 1.0 ) );

    cm.setX0( 0 );
    cm.setY0( 0 );
    cm.setZ0( 0 );
    // x1 is 5, y1 is 6 and z1 is 8

    assert( doubleEq( cm.getDx(), 1.25 ) );
    assert( doubleEq( cm.getDy(), 1.5 ) );
    assert( doubleEq( cm.getDz(), 2.0 ) );

    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testCubeMeshExtendStencil

(

)

Definition at line 726 of file testMesh.cpp.

References CUBOID, MeshCompt::extendStencil(), ChemCompt::getDimensions(), CubeMesh::getMeshDimensions(), CubeMesh::getMeshType(), MeshCompt::getStencilRow(), CubeMesh::innerSetCoords(), and CubeMesh::setPreserveNumEntries().

Referenced by testMesh().

{
    CubeMesh cm0;
    cm0.setPreserveNumEntries( 0 );
    assert( cm0.getMeshType( 0 ) == CUBOID );
    assert( cm0.getMeshDimensions( 0 ) == 3 );
    assert( cm0.getDimensions() == 3 );
    CubeMesh cm1 = cm0;

    vector< double > coords( 9 );
    coords[0] = 0; // X0
    coords[1] = 0; // Y0
    coords[2] = 0; // Z0

    coords[3] = 2; // X1
    coords[4] = 4; // Y1
    coords[5] = 8; // Z1

    coords[6] = 1; // DX
    coords[7] = 1; // DY
    coords[8] = 1; // DZ

    cm0.innerSetCoords( coords );
    coords[2] = 8; // 2x4 face abuts.
    coords[5] = 16;
    cm1.innerSetCoords( coords );

    const double* entry;
    const unsigned int* colIndex;
    unsigned int num = cm0.getStencilRow( 0, &entry, &colIndex );
    assert( num == 3 );
    assert( colIndex[0] == 1 );
    assert( colIndex[1] == 2 );
    assert( colIndex[2] == 8 );

    num = cm0.getStencilRow( 56, &entry, &colIndex );
    assert( num == 3 );
    assert( colIndex[0] == 48 );
    assert( colIndex[1] == 57 );
    assert( colIndex[2] == 58 );

    vector< VoxelJunction > vj;
    for ( unsigned int i = 0; i < 8; ++i ) {
        vj.push_back( VoxelJunction( 56 + i, i ) );
    }
    cm0.extendStencil( &cm1, vj );

    num = cm0.getStencilRow( 56, &entry, &colIndex );
    assert( num == 4 );
    assert( colIndex[0] == 48 );
    assert( colIndex[1] == 57 );
    assert( colIndex[2] == 58 );
    assert( colIndex[3] == 64 );

    for ( unsigned int i = 0; i < 8; ++i ) {
        num = cm0.getStencilRow( 64 + i, &entry, &colIndex );
        assert( num == 1 );
        assert( colIndex[0] == 56 + i );
    }

    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testCubeMeshFillThreeDimSurface

(

)

Definition at line 1447 of file testMesh.cpp.

Referenced by testMesh().

{
    cout << "." << flush;
}

Here is the caller graph for this function:

void testCubeMeshFillTwoDimSurface

(

)

Definition at line 1428 of file testMesh.cpp.

References CubeMesh::innerSetCoords(), CubeMesh::numDims(), CubeMesh::setPreserveNumEntries(), and CubeMesh::surface().

Referenced by testMesh().

{
    CubeMesh cm;
    vector< double > coords( 9, 0.0 );
    coords[3] = 5.0;
    coords[4] = 3.0;
    coords[5] = 1.0;
    coords[6] = coords[7] = coords[8] = 1.0;
    cm.setPreserveNumEntries( false );
    cm.innerSetCoords( coords );
    assert( cm.numDims() == 2 );
    const vector< unsigned int >& surface = cm.surface();
    assert( surface.size() == 15 );
    for ( unsigned int i = 0; i < 15; ++i ) {
        assert( surface[i] == i );
    }
    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testCubeMeshJunctionDiffSizeMesh

(

)

                    14  15
                    12  13

10 11 12 13 14 10 11 8 9 5 6 7 8 9 6 7 4 5 0 1 2 3 4 2 3 0 1

So, junction should be (4,2)(4,4),(9,6),(9,8),(14,10),(14,12)

Definition at line 1560 of file testMesh.cpp.

References CubeMesh::innerSetCoords(), CubeMesh::matchCubeMeshEntries(), CubeMesh::setPreserveNumEntries(), and CubeMesh::surface().

Referenced by testMesh().

{
    CubeMesh cm1;
    vector< double > coords( 9, 0.0 );
    coords[3] = 5.0;
    coords[4] = 3.0;
    coords[5] = 1.0;
    coords[6] = coords[7] = coords[8] = 1.0;
    cm1.setPreserveNumEntries( false );
    cm1.innerSetCoords( coords );
    vector< unsigned int > surface = cm1.surface();
    assert( surface.size() == 15 );

    CubeMesh cm2;
    coords[0] = 5.0;
    coords[1] = -0.5;
    coords[2] = 0.0;
    coords[3] = 7.0;
    coords[4] = 3.5;
    coords[5] = 0.5;
    coords[6] = 1.0;
    coords[7] = 0.5;
    coords[8] = 0.5;
    cm2.setPreserveNumEntries( false );
    cm2.innerSetCoords( coords );
    const vector< unsigned int >& surface2 = cm2.surface();
    assert( surface2.size() == 16 );

    vector< VoxelJunction > ret;
    cm1.matchCubeMeshEntries( &cm2, ret );
    assert( ret.size() == 6 );

    assert( ret[0].first == 4 );
    assert( ret[0].second == 2 );
    assert( ret[1].first == 4 );
    assert( ret[1].second == 4 );
    assert( ret[2].first == 9 );
    assert( ret[2].second == 6 );
    assert( ret[3].first == 9 );
    assert( ret[3].second == 8 );
    assert( ret[4].first == 14 );
    assert( ret[4].second == 10 );
    assert( ret[5].first == 14 );
    assert( ret[5].second == 12 );

    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testCubeMeshJunctionThreeDimSurface

(

)

Definition at line 1620 of file testMesh.cpp.

Referenced by testMesh().

{
    cout << "." << flush;
}

Here is the caller graph for this function:

void testCubeMeshJunctionTwoDimSurface

(

)

                    8   9

10 11 12 13 14 6 7 5 6 7 8 9 4 5 0 1 2 3 4 2 3 0 1

So, junction should be (4,2),(9,4),(14,6)

That was too easy, since the spatial and meshIndices were identical. Now trim the geometries a bit to look like: 6 7 10 11 12 13 - 4 5 5 6 7 8 9 2 3 0 1 2 3 4 - 1

• 0

So, junction should be (9,2) only.

Definition at line 1452 of file testMesh.cpp.

References CubeMesh::getMeshToSpace(), CubeMesh::getSpaceToMesh(), CubeMesh::innerSetCoords(), CubeMesh::matchCubeMeshEntries(), CubeMesh::setMeshToSpace(), CubeMesh::setPreserveNumEntries(), CubeMesh::setSpaceToMesh(), CubeMesh::setSurface(), and CubeMesh::surface().

Referenced by testMesh().

{
    CubeMesh cm1;
    vector< double > coords( 9, 0.0 );
    coords[3] = 5.0;
    coords[4] = 3.0;
    coords[5] = 1.0;
    coords[6] = coords[7] = coords[8] = 1.0;
    cm1.setPreserveNumEntries( false );
    cm1.innerSetCoords( coords );
    vector< unsigned int > surface = cm1.surface();
    assert( surface.size() == 15 );

    CubeMesh cm2;
    coords[0] = 5.0;
    coords[1] = -1.0;
    coords[2] = 0.0;
    coords[3] = 7.0;
    coords[4] = 4.0;
    coords[5] = 1.0;
    coords[6] = coords[7] = coords[8] = 1.0;
    cm2.setPreserveNumEntries( false );
    cm2.innerSetCoords( coords );
    const vector< unsigned int >& surface2 = cm2.surface();
    assert( surface2.size() == 10 );

    vector< VoxelJunction > ret;
    cm1.matchCubeMeshEntries( &cm2, ret );
    assert( ret.size() == 3 );

    assert( ret[0].first == 4 );
    assert( ret[0].second == 2 );
    assert( ret[1].first == 9 );
    assert( ret[1].second == 4 );
    assert( ret[2].first == 14 );
    assert( ret[2].second == 6 );

    // Trimming cm1. At this point we don't assume automatic updates of
    // the m2s, s2m and surface vectors when any of them is changed.
    vector< unsigned int > m2s = cm1.getMeshToSpace();
    assert( m2s.size() == 15 );
    m2s.resize( 14 );
    cm1.setMeshToSpace( m2s );
    vector< unsigned int > s2m = cm1.getSpaceToMesh();
    assert( s2m.size() == 15 );
    s2m[14] = ~0;
    cm1.setSpaceToMesh( s2m );
    surface.resize( 4 );
    // As a shortcut, just assign the places near the junction
    // Note that the indices are spaceIndices.
    surface[0] = 3;
    surface[1] = 4;
    surface[2] = 9;
    surface[3] = 13;
    cm1.setSurface( surface );

    // Trimming cm2.
    m2s = cm2.getMeshToSpace();
    assert( m2s.size() == 10 );
    m2s.resize( 8 );
    m2s[0] = 1;
    for ( unsigned int i = 1; i < 8; ++i )
        m2s[i] = i + 2;
    cm2.setMeshToSpace( m2s );
    s2m.clear();
    s2m.resize( 10, ~0 );
    for ( unsigned int i = 0; i < 8; ++i )
        s2m[ m2s[i] ] = i;
    cm2.setSpaceToMesh( s2m );
    // As a shortcut, just assign the places near the junction
    // Note that the indices are spaceIndices.
    surface[0] = 3;
    surface[1] = 4;
    surface[2] = 6;
    surface[3] = 8;
    cm2.setSurface( surface );

    // Now test it out.
    ret.resize( 0 );
    cm1.matchCubeMeshEntries( &cm2, ret );
    assert( ret.size() == 1 );
    assert( ret[0].first == 9 );
    assert( ret[0].second == 2 );

    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testCubeMeshMultiJunctionTwoD

(

)

The simulated geometry is: D D BBBA CCCCC Here, A is at (0,0,0) to (10,10,10) microns. B is then (-30,0,0) to (0,10,10) microns C is (-30,-10,0) to (20,0,10) microns D is (0,10,0) to (10,30,10) microns.

Definition at line 1636 of file testMesh.cpp.

References CubeMesh::innerSetCoords(), CubeMesh::setPreserveNumEntries(), and CubeMesh::surface().

Referenced by testMesh().

{
    CubeMesh A;
    vector< double > coords( 9, 0.0 );
    coords[3] = 10e-6;
    coords[4] = 10e-6;
    coords[5] = 10e-6;
    coords[6] = coords[7] = coords[8] = 10e-6;
    A.setPreserveNumEntries( false );
    A.innerSetCoords( coords );
    vector< unsigned int > surface = A.surface();
    assert( surface.size() == 1 );
    assert( surface[0] == 0 );

    CubeMesh B;
    coords[0] = -30e-6; coords[1] = 0;      coords[2] = 0;
    coords[3] = 0;      coords[4] = 10e-6;  coords[5] = 10e-6;
    coords[6] =         coords[7] =         coords[8] = 10e-6;
    B.setPreserveNumEntries( false );
    B.innerSetCoords( coords );
    surface = B.surface();
    assert( surface.size() == 3 );
    assert( surface[0] == 0 );
    assert( surface[1] == 1 );
    assert( surface[2] == 2 );

    CubeMesh D;
    coords[0] = 0;      coords[1] = 10e-6;  coords[2] = 0;
    coords[3] = 10e-6;  coords[4] = 30e-6;  coords[5] = 10e-6;
    coords[6] =         coords[7] =         coords[8] = 10e-6;
    D.setPreserveNumEntries( false );
    D.innerSetCoords( coords );
    surface = D.surface();
    assert( surface.size() == 2 );
    assert( surface[0] == 0 );
    assert( surface[1] == 1 );

    CubeMesh C;
    coords[0] = -30e-6; coords[1] = -10e-6; coords[2] = 0;
    coords[3] = 20e-6;  coords[4] = 0;      coords[5] = 10e-6;
    coords[6] =         coords[7] =         coords[8] = 10e-6;
    C.setPreserveNumEntries( false );
    C.innerSetCoords( coords );
    surface = C.surface();
    assert( surface.size() == 5 );
    assert( surface[0] == 0 );
    assert( surface[1] == 1 );
    assert( surface[2] == 2 );
    assert( surface[3] == 3 );
    assert( surface[4] == 4 );


    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testIntersectVoxel

(

)

Here is the geometry of the surface. * is surface, - is empty.

    -***-
    -*---
    -***-


    x***x
    x*32-
    x***x

x is ABUTX
y is ABUTY
z is ABUTZ
2 is 2 points
3 is MULTI

Definition at line 1318 of file testMesh.cpp.

References CubeMesh::ABUTX, checkAbut(), CubeMesh::EMPTY, CubeMesh::MULTI, setIntersectVoxel(), and CubeMesh::SURFACE.

{
    unsigned int nx = 5;
    unsigned int ny = 3;
    unsigned int nz = 1;
    vector< PII > intersect( nx * ny * nz, PII(
                            CubeMesh::EMPTY, CubeMesh::EMPTY ) );
    unsigned int meshIndex = 0;
    setIntersectVoxel( intersect, 1, 0, 0, nx, ny, nz, meshIndex++ );
    setIntersectVoxel( intersect, 2, 0, 0, nx, ny, nz, meshIndex++ );
    setIntersectVoxel( intersect, 3, 0, 0, nx, ny, nz, meshIndex++ );
    setIntersectVoxel( intersect, 1, 1, 0, nx, ny, nz, meshIndex++ );
    setIntersectVoxel( intersect, 1, 2, 0, nx, ny, nz, meshIndex++ );
    setIntersectVoxel( intersect, 2, 2, 0, nx, ny, nz, meshIndex++ );
    setIntersectVoxel( intersect, 3, 2, 0, nx, ny, nz, meshIndex++ );

    assert( intersect[0].first == 0 &&
                    intersect[0].second == CubeMesh::ABUTX );
    assert( intersect[1].first == 0 &&
                    intersect[1].second == CubeMesh::SURFACE );
    assert( intersect[2].first == 1 &&
                    intersect[2].second == CubeMesh::SURFACE );
    assert( intersect[3].first == 2 &&
                    intersect[3].second == CubeMesh::SURFACE );
    assert( intersect[4].first == 2 &&
                    intersect[4].second == CubeMesh::ABUTX );

    assert( intersect[5].first == 3 &&
                    intersect[5].second == CubeMesh::ABUTX );
    assert( intersect[6].first == 3 &&
                    intersect[6].second == CubeMesh::SURFACE );
    assert( intersect[7].first == 1 &&
                    intersect[7].second == CubeMesh::MULTI );
    assert( intersect[8].first == 2 &&
                    intersect[8].second == CubeMesh::MULTI );
    assert( intersect[9].first == EMPTY &&
                    intersect[9].second == CubeMesh::EMPTY );

    assert( intersect[10].first == 4 &&
                    intersect[10].second == CubeMesh::ABUTX );
    assert( intersect[11].first == 4 &&
                    intersect[11].second == CubeMesh::SURFACE );
    assert( intersect[12].first == 5 &&
                    intersect[12].second == CubeMesh::SURFACE );
    assert( intersect[13].first == 6 &&
                    intersect[13].second == CubeMesh::SURFACE );
    assert( intersect[14].first == 6 &&
                    intersect[14].second == CubeMesh::ABUTX );

    // Next: test out checkAbut.
    vector< VoxelJunction > ret;
    checkAbut( intersect, 0, 0, 0, nx, ny, nz, 1234, ret );
    assert( ret.size() == 1 );
    assert( ret[0].first == 0 && ret[0].second == 1234 );
    ret.clear();
    // The ones below are either SURFACE or EMPTY and should not add points
    checkAbut( intersect, 1, 0, 0, nx, ny, nz, 1234, ret );
    checkAbut( intersect, 2, 0, 0, nx, ny, nz, 1234, ret );
    checkAbut( intersect, 3, 0, 0, nx, ny, nz, 1234, ret );
    checkAbut( intersect, 1, 1, 0, nx, ny, nz, 1234, ret );
    checkAbut( intersect, 4, 1, 0, nx, ny, nz, 1234, ret );
    checkAbut( intersect, 1, 2, 0, nx, ny, nz, 1234, ret );
    checkAbut( intersect, 2, 2, 0, nx, ny, nz, 1234, ret );
    checkAbut( intersect, 3, 2, 0, nx, ny, nz, 1234, ret );
    assert( ret.size() == 0 );
    checkAbut( intersect, 2, 1, 0, nx, ny, nz, 9999, ret );
    assert( ret.size() == 3 );
    assert( ret[0].first == 3 && ret[0].second == 9999 );
    assert( ret[1].first == 1 && ret[1].second == 9999 );
    assert( ret[2].first == 5 && ret[1].second == 9999 );
    ret.clear();
    checkAbut( intersect, 3, 1, 0, nx, ny, nz, 8888, ret );
    assert( ret.size() == 2 );
    assert( ret[0].first == 2 && ret[0].second == 8888 );
    assert( ret[1].first == 6 && ret[1].second == 8888 );
    ret.clear();
    checkAbut( intersect, 4, 0, 0, nx, ny, nz, 7777, ret );
    checkAbut( intersect, 0, 1, 0, nx, ny, nz, 6666, ret );
    checkAbut( intersect, 0, 2, 0, nx, ny, nz, 5555, ret );
    checkAbut( intersect, 4, 2, 0, nx, ny, nz, 4444, ret );
    assert( ret.size() == 4 );
    assert( ret[0].first == 2 && ret[0].second == 7777 );
    assert( ret[1].first == 3 && ret[1].second == 6666 );
    assert( ret[2].first == 4 && ret[2].second == 5555 );
    assert( ret[3].first == 6 && ret[3].second == 4444 );

    cout << "." << flush;
}

Here is the call graph for this function:

void testMesh

(

)

Definition at line 2138 of file testMesh.cpp.

References testCubeMesh(), testCubeMeshExtendStencil(), testCubeMeshFillThreeDimSurface(), testCubeMeshFillTwoDimSurface(), testCubeMeshJunctionDiffSizeMesh(), testCubeMeshJunctionThreeDimSurface(), testCubeMeshJunctionTwoDimSurface(), testCubeMeshMultiJunctionTwoD(), testVec(), and testVolScaling().

Referenced by nonMpiTests().

{
    testVec();
    testVolScaling();
    // testCylBase();
    // testNeuroNode();
    // testNeuroNodeTree();
    // testCylMesh();
    // testMidLevelCylMesh();
    testCubeMesh();
    testCubeMeshExtendStencil();
    // testReMesh(); // Waiting to have pool subdivision propagate.
    // testNeuroMeshLinear();
    // testNeuroMeshBranching();
    // testIntersectVoxel();
    testCubeMeshFillTwoDimSurface();
    testCubeMeshFillThreeDimSurface();
    testCubeMeshJunctionTwoDimSurface();
    testCubeMeshJunctionThreeDimSurface();
    testCubeMeshJunctionDiffSizeMesh();
    testCubeMeshMultiJunctionTwoD();
    // testSpineEntry();
    // testSpineAndPsdMesh();
    // testCellPortion();
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testReMesh

(

)

Tests scaling of volume and # of entries upon mesh resize

Definition at line 792 of file testMesh.cpp.

References Eref::data(), Shell::doCreate(), Shell::doDelete(), doubleEq(), Id::eref(), Field< A >::get(), NA, Field< A >::set(), and SetGet2< A1, A2 >::set().

{
    Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
    Id base = s->doCreate( "Neutral", Id(), "base", 1 );

    Id cube = s->doCreate( "CubeMesh", base, "cube", 1 );
    bool ret = SetGet2< double, unsigned int >::set(
        cube, "buildDefaultMesh", 1.0, 1 );
    assert( ret );
    unsigned int vol = Field< double >::get( cube, "volume" );
    assert( doubleEq( vol, 1.0 ) );
    Id pool = s->doCreate( "Pool", cube, "pool", 1 );
    Id mesh( "/base/cube/mesh" );
    assert( mesh != Id() );

    // 1 millimolar in 1 m^3 is 1 mole per liter.
    double linsize = Field< double >::get( pool, "volume" );
    assert( doubleEq( linsize, 1.0 ) );

    ret = Field< double >::set( pool, "conc", 1 );
    assert( ret );
    double n = Field< double >::get( pool, "n" );
    assert( doubleEq( n, NA ) );

    ret = SetGet2< double, unsigned int >::set(
        cube, "buildDefaultMesh", 1.0e-3, 1 );
    Field< double >::set( pool, "conc", 1 );
    n = Field< double >::get( pool, "n" );
    assert( doubleEq( n, NA / 1000.0 ) );

    // Next we do the remeshing.
    double x = 1.234;
    Field< double >::set( pool, "concInit", x );
    ret = SetGet2< double, unsigned int >::set(
        cube, "buildDefaultMesh", 1, 8 );
    // This is nasty, needs the change to propagate through messages.
    linsize = Field< double >::get( pool, "volume" );
    assert( doubleEq( linsize, 0.125 ) );

    n = Field< double >::get( ObjId( pool, 0 ), "concInit" );
    assert( doubleEq( n, x ) );
    n = Field< double >::get( ObjId( pool, 7 ), "concInit" );
    assert( doubleEq( n, x ) );
    n = Field< double >::get( ObjId( pool, 0 ), "nInit" );
    assert( doubleEq( n, x * NA / 8.0 ) );
    n = Field< double >::get( ObjId( pool, 7 ), "nInit" );
    assert( doubleEq( n, x * NA / 8.0 ) );
    n = Field< double >::get( ObjId( pool, 0 ), "conc" );
    assert( doubleEq( n, x ) );
    n = Field< double >::get( ObjId( pool, 7 ), "conc" );
    assert( doubleEq( n, x ) );

    s->doDelete( base );
    cout << "." << flush;
}

Here is the call graph for this function:

void testVec

(

)

Definition at line 1691 of file testMesh.cpp.

References Vec::crossProduct(), Vec::dotProduct(), doubleEq(), and Vec::orthogonalAxes().

Referenced by testMesh().

{
    Vec i( 1, 0, 0 );
    Vec j( 0, 1, 0 );
    Vec k( 0, 0, 1 );

    assert( doubleEq( i.dotProduct( j ), 0.0 ) );
    assert( doubleEq( j.dotProduct( k ), 0.0 ) );
    assert( doubleEq( j.dotProduct( j ), 1.0 ) );

    assert( i.crossProduct( j ) == k );

    Vec u, v;
    i.orthogonalAxes( u, v );
    assert( u == j );
    assert( v == k );

    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void testVolScaling

(

)

This tests how volume changes in a mesh propagate to all child pools, reacs, and enzymes.

Definition at line 34 of file testMesh.cpp.

References Eref::data(), Shell::doCreate(), Shell::doDelete(), doubleEq(), Id::eref(), Field< A >::get(), makeReacTest(), and Field< A >::set().

Referenced by testMesh().

{
    extern Id makeReacTest();
    Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );

    Id kin = makeReacTest();
    Id subCompt =  s->doCreate( "CubeMesh", kin, "subCompt", 1 );
    Field< double >::set( subCompt, "volume", 1e-16 );
    Id SP = s->doCreate( "Pool", subCompt, "SP", 1 );
    Field< double >::set( SP, "concInit", 2 );

    vector< double > n;
    n.push_back( Field< double >::get( ObjId( "/kinetics/A" ), "nInit" ) );
    n.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool" ), "nInit" ) );
    n.push_back( Field< double >::get( ObjId( "/kinetics/r1" ), "numKf" ));
    n.push_back( Field< double >::get( ObjId( "/kinetics/r1" ), "numKb" ));
    n.push_back( Field< double >::get( ObjId( "/kinetics/r2" ), "numKf" ));
    n.push_back( Field< double >::get( ObjId( "/kinetics/r2" ), "numKb" ));
    n.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k1" ) );
    n.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k2" ) );
    n.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k3" ) );
    n.push_back( Field< double >::get( ObjId( "/kinetics/e2Pool/e2" ), "Km" ) );
    n.push_back( Field< double >::get( ObjId( "/kinetics/e2Pool/e2" ), "kcat" ) );
    n.push_back( Field< double >::get( SP, "nInit" ) );

    double vol = Field< double >::get( kin, "volume" );

    vol *= 10;
    Field< double >::set( kin, "volume", vol );

    vector< double > m;
    m.push_back( Field< double >::get( ObjId( "/kinetics/A" ), "nInit" ) );
    m.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool" ), "nInit" ) );
    m.push_back( Field< double >::get( ObjId( "/kinetics/r1" ), "numKf" ));
    m.push_back( Field< double >::get( ObjId( "/kinetics/r1" ), "numKb" ));
    m.push_back( Field< double >::get( ObjId( "/kinetics/r2" ), "numKf" ));
    m.push_back( Field< double >::get( ObjId( "/kinetics/r2" ), "numKb" ));
    m.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k1" ) );
    m.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k2" ) );
    m.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k3" ) );
    m.push_back( Field< double >::get( ObjId( "/kinetics/e2Pool/e2" ), "Km" ) );
    m.push_back( Field< double >::get( ObjId( "/kinetics/e2Pool/e2" ), "kcat" ) );
    m.push_back( Field< double >::get( SP, "nInit" ) );

    assert( doubleEq( n[0] * 10, m[0] ) );  // A nInit
    assert( doubleEq( n[1] * 10, m[1] ) );  // e1Pool
    assert( doubleEq( n[2] / 10, m[2] ) );  // r1 numKf
    assert( doubleEq( n[3], m[3] ) );   // r1 numKb
    assert( doubleEq( n[4] / 10, m[4] ) );  // r2 numKf
    assert( doubleEq( n[5], m[5] ) );   // r2 numKb
    assert( doubleEq( n[6] / 10, m[6] ) );  // e1 k1
    assert( doubleEq( n[7], m[7] ) );   // e1 k2
    assert( doubleEq( n[8], m[8] ) );   // e1 k3
    assert( doubleEq( n[9], m[9] ) );   // e2 Km
    assert( doubleEq( n[10], m[10] ) ); // e2 kcat
    assert( doubleEq( n[11], m[11] ) ); // SP nInit

    s->doDelete( kin );
    cout << "." << flush;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

const unsigned int ABUT = ~2

static

Definition at line 1302 of file testMesh.cpp.

const unsigned int EMPTY = ~0

static

Definition at line 1300 of file testMesh.cpp.

const unsigned int MULTI = ~3

static

Definition at line 1303 of file testMesh.cpp.

const unsigned int SURFACE = ~1

static

Definition at line 1301 of file testMesh.cpp.