C++ Shell::doCreate方法代码示例

void testRunKsolve()
{
	double simDt = 0.1;
	// double plotDt = 0.1;
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	Id kin = makeReacTest();
	Id ksolve = s->doCreate( "Ksolve", kin, "ksolve", 1 );
	Id stoich = s->doCreate( "Stoich", ksolve, "stoich", 1 );
	Field< Id >::set( stoich, "compartment", kin );
	Field< Id >::set( stoich, "ksolve", ksolve );
	Field< string >::set( stoich, "path", "/kinetics/##" );
	s->doUseClock( "/kinetics/ksolve", "process", 4 ); 
	s->doSetClock( 4, simDt );

	s->doReinit();
	s->doStart( 20.0 );
	Id plots( "/kinetics/plots" );
	for ( unsigned int i = 0; i < 7; ++i ) {
		stringstream ss;
		ss << "plot." << i;
		SetGet2< string, string >::set( ObjId( plots, i ), "xplot", 
						"tsr2.plot", ss.str() );
	}
	s->doDelete( kin );
	cout << "." << flush;
}

void testReacVolumeScaling()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	Id comptId = shell->doCreate( "CubeMesh", Id(), "cube", 1 );
	Id meshId( comptId.value() + 1 );
	Id subId = shell->doCreate( "Pool", comptId, "sub", 1 );
	Id prdId = shell->doCreate( "Pool", comptId, "prd", 1 );
	Id reacId = shell->doCreate( "Reac", comptId, "reac", 1 );

	double vol1 = 1e-15;

	ObjId mid = shell->doAddMsg( "OneToOne", 
		subId, "requestVolume", meshId, "get_volume" );
	assert( mid != ObjId() );
	mid = shell->doAddMsg( "OneToOne", 
		prdId, "requestVolume", meshId, "get_volume" );
	assert( mid != ObjId() );

	vector< double > coords( 9, 10.0e-6 );
	coords[0] = coords[1] = coords[2] = 0;

	Field< vector< double > >::set( comptId, "coords", coords );

	double volume = Field< double >::get( comptId, "volume" );
	assert( doubleEq( volume, vol1 ) );

	ObjId ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );
	ret = shell->doAddMsg( "Single", reacId, "prd", prdId, "reac" );
	assert( ret != ObjId() );

	Field< double >::set( reacId, "Kf", 2 );
	Field< double >::set( reacId, "Kb", 3 );
	double x = Field< double >::get( reacId, "kf" );
	assert( doubleEq( x, 2 ) );
	x = Field< double >::get( reacId, "kb" );
	assert( doubleEq( x, 3 ) );
	
	ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );
	double conv = 1.0 / ( NA * vol1 );
	x = Field< double >::get( reacId, "kf" );
	assert( doubleEq( x, 2 * conv ) );
	x = Field< double >::get( reacId, "kb" );
	assert( doubleEq( x, 3 ) );

	ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );
	ret = shell->doAddMsg( "Single", reacId, "prd", prdId, "reac" );
	assert( ret != ObjId() );
	x = Field< double >::get( reacId, "kf" );
	assert( doubleEq( x, 2 * conv * conv ) );
	x = Field< double >::get( reacId, "kb" );
	assert( doubleEq( x, 3 * conv ) );

	shell->doDelete( comptId );
	cout << "." << flush;
}

void testMMenzProcess()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	//////////////////////////////////////////////////////////////////////
	// This set is the test kinetic calculation using MathFunc
	//////////////////////////////////////////////////////////////////////
	Id nid = shell->doCreate( "Neutral", Id(), "n", 1 );
	//////////////////////////////////////////////////////////////////////
	// This set is the reference kinetic calculation using MMEnz
	//////////////////////////////////////////////////////////////////////
	Id pid = shell->doCreate( "Pool", nid, "p", 1 ); // substrate
	Id qid = shell->doCreate( "Pool", nid, "q", 1 );	// enz mol
	Id rid = shell->doCreate( "Pool", nid, "r", 1 ); // product
	Id mmid = shell->doCreate( "MMenz", nid, "mm", 1 ); // mmenz

	Id tabid2 = shell->doCreate( "Table", nid, "tab2", 1 ); //output plot

	Field< double >::set( mmid, "Km", 1.0 );
	Field< double >::set( mmid, "kcat", 1.0 );
	Field< double >::set( pid, "nInit", 1.0 );
	Field< double >::set( qid, "nInit", 1.0 );
	Field< double >::set( rid, "nInit", 0.0 );

	shell->doAddMsg( "Single", ObjId( mmid ), "sub", ObjId( pid ), "reac" );
	shell->doAddMsg( "Single", ObjId( mmid ), "prd", ObjId( rid ), "reac" );
	shell->doAddMsg( "Single", ObjId( qid ), "nOut", ObjId( mmid ), "enzDest" );
	shell->doAddMsg( "Single", ObjId( pid ), "nOut", ObjId( tabid2 ), "input" );
	shell->doSetClock( 0, 0.01 );
	shell->doSetClock( 1, 0.01 );
	shell->doUseClock( "/n/mm,/n/tab2", "process", 0 );
	shell->doUseClock( "/n/#[ISA=Pool]", "process", 1 );
	
	//////////////////////////////////////////////////////////////////////
	// Now run models and compare outputs
	//////////////////////////////////////////////////////////////////////

	shell->doReinit();
	shell->doStart( 10 );

	vector< double > vec = Field< vector< double > >::get( tabid2, "vec" );
	assert( vec.size() == 1001 );
	for ( unsigned int i = 0; i < vec.size(); ++i ) {
		double t = 0.01 * i;
		double et = estT( vec[i] );
		assert( doubleApprox( t, et ) );
	}

	shell->doDelete( nid );
	cout << "." << flush;
}

void testTaperingCylDiffn()
{
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = 25e-6;
	double r0 = 2e-6;
	double r1 = 1e-6;
	double diffLength = 1e-6; // 1e-6 is the highest dx for which error is OK
	double runtime = 10.0;
	double dt = 0.1; // 0.2 is the highest dt for which the error is in bounds
	double diffConst = 1.0e-12; 
	// Should set explicitly, currently during creation of DiffPoolVec
	//double diffConst = 1.0e-12; 
	Id model = s->doCreate( "Neutral", Id(), "model", 1 );
	Id cyl = s->doCreate( "CylMesh", model, "cyl", 1 );
	Field< double >::set( cyl, "r0", r0 );
	Field< double >::set( cyl, "r1", r1 );
	Field< double >::set( cyl, "x0", 0 );
	Field< double >::set( cyl, "x1", len );
	Field< double >::set( cyl, "diffLength", diffLength );
	unsigned int ndc = Field< unsigned int >::get( cyl, "numMesh" );
	assert( ndc == static_cast< unsigned int >( round( len / diffLength )));
	Id pool = s->doCreate( "Pool", cyl, "pool", 1 );
	Field< double >::set( pool, "diffConst", diffConst );

	Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
	Field< Id >::set( dsolve, "compartment", cyl );
	s->doUseClock( "/model/dsolve", "process", 1 );
	s->doSetClock( 1, dt );
	// Next: build by setting the path of the dsolve.
	Field< string >::set( dsolve, "path", "/model/cyl/pool" );
	// Then find a way to test it.
	assert( pool.element()->numData() == ndc );
	Field< double >::set( ObjId( pool, 0 ), "nInit", 1.0 );

	s->doReinit();
	s->doStart( runtime );

	double myTot = 0.0;
	vector< double > poolVec;
   	Field< double >::getVec( pool, "n", poolVec );
	for ( unsigned int i = 0; i < poolVec.size(); ++i ) {
		myTot += poolVec[i];
	} 
	assert( doubleEq( myTot, 1.0 ) );

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

/// crate test object and push it into the container vector
Id create_testobject(vector<Id> & container, string classname, Id parent, string name, vector<int> dims)
{
    Shell * shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() );
    Id nid = shell->doCreate( classname, parent, name, dims );
    container.push_back(nid);
    return nid;
}

/* Enzymatic Reaction  */
void SbmlReader::setupEnzymaticReaction( const EnzymeInfo & einfo,string enzname, const map< string, Id > &molSidcmptMIdMap,string name1) {
    string enzPool = einfo.enzyme;

    Id comptRef = molSidcmptMIdMap.find(enzPool)->second; //gives compartment of sp
    Id meshEntry = Neutral::child( comptRef.eref(), "mesh" );
    Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );

    //Creating enz pool to enzyme site
    Id enzPoolId = molSidMIdMap_.find(enzPool)->second;
    Id enzyme_ = shell->doCreate("Enz", enzPoolId, name1, 1);
    //shell->doAddMsg( "Single", meshEntry, "remeshReacs", enzyme_, "remesh");

    Id complex = einfo.complex;
    //Moving enzyme site under enzyme
    shell->doMove(complex,enzyme_);
    shell->doAddMsg("OneToAll",enzyme_,"cplx",complex,"reac");

    shell->doAddMsg("OneToOne",enzyme_,"enz",enzPoolId,"reac");

    vector< Id >::const_iterator sub_itr;
    for ( sub_itr = einfo.substrates.begin(); sub_itr != einfo.substrates.end(); sub_itr++ ) {
        Id S = (*sub_itr);
        Id b = shell->doAddMsg( "OneToOne", enzyme_, "sub" ,S , "reac" );
    }

    vector< Id >::const_iterator prd_itr;
    for ( prd_itr = einfo.products.begin(); prd_itr != einfo.products.end(); prd_itr++ ) {
        Id P = (*prd_itr);
        shell->doAddMsg ("OneToOne",enzyme_,"prd", P,"reac");
    }
    // populate k3,k2,k1 in this order only.
    Field < double > :: set( enzyme_, "k3", einfo.k3 );
    Field < double > :: set( enzyme_, "k2", einfo.k2 );
    Field < double > :: set( enzyme_, "k1", einfo.k1 );
}

void testPoolVolumeScaling()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	Id comptId = shell->doCreate( "CylMesh", Id(), "cyl", 1 );
	Id meshId( comptId.value() + 1 );
	Id poolId = shell->doCreate( "Pool", comptId, "pool", 1 );

	ObjId mid = shell->doAddMsg( "OneToOne", 
		ObjId( poolId, 0 ), "requestVolume",
		ObjId( meshId, 0 ), "get_volume" );

	assert( mid != ObjId() );

	vector< double > coords( 9, 0.0 );
	double x1 = 100e-6;
	double r0 = 10e-6;
	double r1 = 5e-6;
	double lambda = x1;
	coords[3] = x1;
	coords[6] = r0;
	coords[7] = r1;
	coords[8] = lambda;

	Field< vector< double > >::set( comptId, "coords", coords );

	double volume = Field< double >::get( poolId, "volume" );
	assert( doubleEq( volume, PI * x1 * (r0+r1) * (r0+r1) / 4.0 ) );

	Field< double >::set( poolId, "n", 400 );
	double volscale = 1 / ( NA * volume );
	double conc = Field< double >::get( poolId, "conc" );
	assert( doubleEq( conc, 400 * volscale ) );
	Field< double >::set( poolId, "conc", 500 * volscale );
	double n = Field< double >::get( poolId, "n" );
	assert( doubleEq( n, 500 ) );

	Field< double >::set( poolId, "nInit", 650 );
	double concInit = Field< double >::get( poolId, "concInit" );
	assert( doubleEq( concInit, 650 * volscale ) );
	Field< double >::set( poolId, "concInit", 10 * volscale );
	n = Field< double >::get( poolId, "nInit" );
	assert( doubleEq( n, 10 ) );

	shell->doDelete( comptId );
	cout << "." << flush;
}

void SbmlReader::getRules() {
    unsigned int nr = model_->getNumRules();
    //if (nr > 0)
    //  cout << "\n ##### Need to populate funcpool and sumtotal which is pending due to equations \n";
    Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
    for ( unsigned int r = 0; r < nr; r++ ) {

        Rule * rule = model_->getRule(r);
        bool assignRule = rule->isAssignment();
        if ( assignRule ) {
            string rule_variable = rule->getVariable();
            map< string,Id >::iterator v_iter;
            map< string,Id >::iterator m_iter;
            v_iter = molSidMIdMap_.find( rule_variable );
            if (v_iter != molSidMIdMap_.end()) {
                Id rVariable = molSidMIdMap_.find(rule_variable)->second;
                string rstring =molSidMIdMap_.find(rule_variable)->first;
                Id sumId = shell->doCreate( "SumFunc", rVariable, "func", 1 );
                rVariable.element()->zombieSwap( FuncPool::initCinfo() );
                ObjId ret = shell->doAddMsg( "single",
                                             ObjId( sumId, 0 ), "output",
                                             ObjId( rVariable, 0 ), "input" );
                assert( ret != ObjId() );
                const ASTNode * ast = rule->getMath();
                vector< string > ruleMembers;
                ruleMembers.clear();
                printMembers( ast,ruleMembers );
                for ( unsigned int rm = 0; rm < ruleMembers.size(); rm++ ) {
                    m_iter = molSidMIdMap_.find( ruleMembers[rm] );
                    if ( m_iter != molSidMIdMap_.end() ) {
                        Id rMember = molSidMIdMap_.find(ruleMembers[rm])->second;
                        ObjId ret = shell->doAddMsg( "single",
                                                     ObjId( rMember, 0 ), "nOut",
                                                     ObjId( sumId, 0 ), "input" );
                        string test = molSidMIdMap_.find(ruleMembers[rm])->first;
                    } else {
                        cerr << "SbmlReader::getRules: Assignment rule member is not a species" << endl;
                        // In assignment rule there are constants instead of molecule which is yet to deal in moose.
                        errorFlag_ = true;
                    }
                }
            }
        }
        bool rateRule = rule->isRate();
        if ( rateRule ) {
            cout << "warning : for now Rate Rule is not handled " << endl;
            errorFlag_ = true;
        }
        bool  algebRule = rule->isAlgebraic();
        if ( algebRule ) {
            cout << "warning: for now Algebraic Rule is not handled" << endl;
            errorFlag_ = true;
        }
    }
}

void testRunGsolve()
{
	double simDt = 0.1;
	// double plotDt = 0.1;
	Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
	Id kin = makeReacTest();
	double volume = 1e-21;
	Field< double >::set( kin, "volume", volume );
	Field< double >::set( ObjId( "/kinetics/A" ), "concInit", 2 );
	Field< double >::set( ObjId( "/kinetics/e1Pool" ), "concInit", 1 );
	Field< double >::set( ObjId( "/kinetics/e2Pool" ), "concInit", 1 );
	Id e1( "/kinetics/e1Pool/e1" );
	Field< double >::set( e1, "Km", 5 );
	Field< double >::set( e1, "kcat", 1 );
	vector< double > stim( 100, 0.0 );
	for ( unsigned int i = 0; i< 100; ++i ) {
		stim[i] = volume * NA * (1.0 + sin( i * 2.0 * PI / 100.0 ) );
	}
	Field< vector< double > >::set( ObjId( "/kinetics/tab" ), "vector", stim );


	Id gsolve = s->doCreate( "Gsolve", kin, "gsolve", 1 );
	Id stoich = s->doCreate( "Stoich", gsolve, "stoich", 1 );
	Field< Id >::set( stoich, "compartment", kin );
	Field< Id >::set( stoich, "ksolve", gsolve );
	
	Field< string >::set( stoich, "path", "/kinetics/##" );
	s->doUseClock( "/kinetics/gsolve", "process", 4 ); 
	s->doSetClock( 4, simDt );

	s->doReinit();
	s->doStart( 20.0 );
	Id plots( "/kinetics/plots" );
	for ( unsigned int i = 0; i < 7; ++i ) {
		stringstream ss;
		ss << "plot." << i;
		SetGet2< string, string >::set( ObjId( plots, i ), "xplot", 
						"tsr3.plot", ss.str() );
	}
	s->doDelete( kin );
	cout << "." << flush;
}

// See what Element::getNeighbors does with 2 sub <----> prd.
void testTwoReacGetNeighbors()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	Id comptId = shell->doCreate( "CubeMesh", Id(), "cube", 1 );
	Id meshId( comptId.value() + 1 );
	Id subId = shell->doCreate( "Pool", comptId, "sub", 1 );
	Id prdId = shell->doCreate( "Pool", comptId, "prd", 1 );
	Id reacId = shell->doCreate( "Reac", comptId, "reac", 1 );

	ObjId mid = shell->doAddMsg( "OneToOne", 
		subId, "requestVolume", meshId, "get_volume" );
	assert( mid != ObjId() );
	mid = shell->doAddMsg( "OneToOne", 
		prdId, "requestVolume", meshId, "get_volume" );
	assert( mid != ObjId() );

	ObjId ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );
	ret = shell->doAddMsg( "Single", reacId, "sub", subId, "reac" );
	assert( ret != ObjId() );

	ret = shell->doAddMsg( "Single", reacId, "prd", prdId, "reac" );
	assert( ret != ObjId() );

	vector< Id > pools;
	unsigned int num = reacId.element()->getNeighbors( pools, 
		Reac::initCinfo()->findFinfo( "toSub" ) );
	assert( num == 2 );
	assert( pools[0] == subId );
	assert( pools[1] == subId );

	pools.clear();
	num = reacId.element()->getNeighbors( pools, 
		Reac::initCinfo()->findFinfo( "sub" ) );
	assert( num == 2 );
	assert( pools[0] == subId );
	assert( pools[1] == subId );

	shell->doDelete( comptId );
	cout << "." << flush;
}

//////////////////////////////////////////////////////////////////
// The read functions.
//////////////////////////////////////////////////////////////////
Id  makeStandardElements( Id pa, const string& modelname )
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	//cout << " kkit read " << pa << " " << modelname << " "<< MooseGlobal;
	string modelPath = pa.path() + "/" + modelname;
	if ( pa == Id() )
		modelPath = "/" + modelname;
	Id mgr( modelPath );
	if ( mgr == Id() )
		mgr = shell->doCreate( "Neutral", pa, modelname, 1, MooseGlobal );
	Id kinetics( modelPath + "/kinetics" );
	if ( kinetics == Id() ) {
		kinetics = 
		shell->doCreate( "CubeMesh", mgr, "kinetics", 1,  MooseGlobal );
		SetGet2< double, unsigned int >::set( kinetics, "buildDefaultMesh", 1e-15, 1 );
	}
	assert( kinetics != Id() );

	Id graphs = shell->doCreate( "Neutral", mgr, "graphs", 1, MooseGlobal);
	assert( graphs != Id() );
	Id moregraphs = shell->doCreate( "Neutral", mgr, "moregraphs", 1, MooseGlobal );

	Id geometry = shell->doCreate( "Neutral", mgr, "geometry", 1, MooseGlobal );
	assert( geometry != Id() );

	Id groups = 
		shell->doCreate( "Neutral", mgr, "groups", 1, MooseGlobal );
	assert( groups != Id() );
	return mgr;
}

void testSpikeGen()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	Id sid = shell->doCreate( "SpikeGen", Id(), "spike", 1, MooseGlobal );
	SpikeGen& sg = *( reinterpret_cast< SpikeGen* >( sid.eref().data() ) );

	Eref er( sid.eref() );
	ProcInfo p;
	p.dt = 0.001;
	p.currTime = 0.0;
	sg.setThreshold( 1.0 );
	sg.setRefractT( 0.005 );

	sg.reinit( er, &p );
	sg.handleVm( 0.5 );
	sg.process( er, &p );
	assert( !sg.getFired() );
	p.currTime += p.dt;

	sg.handleVm( 0.999 );
	sg.process( er, &p );
	assert( !sg.getFired() );
	p.currTime += p.dt;

	sg.handleVm( 1.001 );
	sg.process( er, &p );
	assert( sg.getFired() );
	p.currTime += p.dt;

	sg.handleVm( 0.999 );
	sg.process( er, &p );
	assert( !sg.getFired() );
	p.currTime += p.dt;

	sg.handleVm( 2.0 ); // Too soon, refractory
	sg.process( er, &p );
	assert( !sg.getFired() );

	p.currTime += 0.005; // Now post-refractory
	sg.handleVm( 2.0 ); // Now not refractory
	sg.process( er, &p );
	assert( sg.getFired() );

	sid.destroy();
	cout << "." << flush;
}

/*    set up Michalies Menten reaction  */
void SbmlReader::setupMMEnzymeReaction( Reaction * reac,string rid,string rname,const map< string, Id > &molSidcmptMIdMap ) {
    string::size_type loc = rid.find( "_MM_Reaction_" );
    if( loc != string::npos ) {
        int strlen = rid.length();
        rid.erase( loc,strlen-loc );
    }
    unsigned int num_modifr = reac->getNumModifiers();
    for ( unsigned int m = 0; m < num_modifr; m++ ) {
        const ModifierSpeciesReference* modifr=reac->getModifier( m );
        std::string sp = modifr->getSpecies();
        Id enzyme_;
        Id E = molSidMIdMap_.find(sp)->second;

        Id comptRef = molSidcmptMIdMap.find(sp)->second; //gives compartment of sp
        Id meshEntry = Neutral::child( comptRef.eref(), "mesh" );
        Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
        enzyme_ = shell->doCreate("MMenz",E,rname,1);
        //shell->doAddMsg( "Single", meshEntry, "remeshReacs", enzyme_, "remesh");
        shell->doAddMsg("Single",E,"nOut",enzyme_,"enzDest");

        KineticLaw * klaw=reac->getKineticLaw();
        vector< double > rate;
        rate.clear();
        getKLaw( klaw,true,rate );
        if ( errorFlag_ )
            return;
        else if ( !errorFlag_ ) {
            for ( unsigned int rt = 0; rt < reac->getNumReactants(); rt++ ) {
                const SpeciesReference* rct = reac->getReactant( rt );
                sp=rct->getSpecies();
                Id S = molSidMIdMap_.find(sp)->second;
                shell->doAddMsg( "OneToOne", enzyme_, "sub" ,S , "reac" );
            }
            for ( unsigned int pt = 0; pt < reac->getNumProducts(); pt++ ) {
                const SpeciesReference* pdt = reac->getProduct(pt);
                sp = pdt->getSpecies();
                Id P = molSidMIdMap_.find(sp)->second;
                shell->doAddMsg( "OneToOne", enzyme_, "prd" ,P, "reac" );
            }
            Field < double > :: set( enzyme_, "kcat", rate[0] );
            Field < double > :: set( enzyme_, "numKm", rate[1] );
        }

    }
}

static Id makeCompt( Id parent, 
		const SwcSegment& seg, const SwcSegment& pa,
		double RM, double RA, double CM,
		unsigned int i, unsigned int j	)
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	double len = seg.radius() * 2.0;
	string name = "soma";
	Id compt;
	double x0, y0, z0;
	if ( seg.parent() != ~0U ) {
		len = seg.distance( pa );
		stringstream ss;
		ss << SwcSegment::typeName[ seg.type() ] << "_" << i << "_" << j;
		name = ss.str();
		x0 = pa.vec().a0();
		y0 = pa.vec().a1();
		z0 = pa.vec().a2();
	} else {
		x0 = seg.vec().a0() - len;
		y0 = seg.vec().a1();
		z0 = seg.vec().a2();
	}
	assert( len > 0.0 );
	compt = shell->doCreate( "Compartment", parent, name, 1 );
	Eref er = compt.eref();
	moose::CompartmentBase *cptr = reinterpret_cast< moose::CompartmentBase* >(
					compt.eref().data() );
	double xa = seg.radius() * seg.radius() * PI * 1e-12;
	len *= 1e-6;
	double dia = seg.radius() * 2.0e-6;
	cptr->setRm( er, RM / ( len * dia * PI ) );
	cptr->setRa( er, RA * len / xa );
	cptr->setCm( er, CM * ( len * dia * PI ) );
	cptr->setDiameter( dia );
	cptr->setLength( len );
	cptr->setX0( x0 * 1e-6 );
	cptr->setY0( y0 * 1e-6 );
	cptr->setZ0( z0 * 1e-6 );
	cptr->setX( seg.vec().a0() * 1e-6 );
	cptr->setY( seg.vec().a1() * 1e-6 );
	cptr->setZ( seg.vec().a2() * 1e-6 );
	return compt;
}

void testMMenz()
{
	Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
	Id mmid = shell->doCreate( "MMenz", Id(), "mm", 1 ); // mmenz
	MMenz m;
	ProcInfo p;

	m.vSetKm( mmid.eref(), 5.0 );
	m.vSetKcat( mmid.eref(), 4.0 );
	m.vReinit( mmid.eref(), &p );
	m.vSub( 2 );
	m.vEnz( 3 );
	assert( doubleEq( m.vGetKm( mmid.eref() ), 5.0 ) );
	assert( doubleEq( m.vGetKcat( mmid.eref() ), 4.0 ) );
	m.vProcess( mmid.eref(), &p );

	shell->doDelete( mmid );
	cout << "." << flush;
}