MOOSE - Multiscale Object Oriented Simulation Environment
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
testMesh.cpp
Go to the documentation of this file.
1 /**********************************************************************
2 ** This program is part of 'MOOSE', the
3 ** Messaging Object Oriented Simulation Environment.
4 ** Copyright (C) 2011 Upinder S. Bhalla. and NCBS
5 ** It is made available under the terms of the
6 ** GNU Lesser General Public License version 2.1
7 ** See the file COPYING.LIB for the full notice.
8 **********************************************************************/
9 
10 #include "header.h"
11 #include "SparseMatrix.h"
12 #include "../shell/Shell.h"
13 #include "Boundary.h"
14 #include "MeshEntry.h"
15 // #include "Stencil.h"
16 #include "ChemCompt.h"
17 #include "MeshCompt.h"
18 #include "CubeMesh.h"
19 #include "CylBase.h"
20 #include "NeuroNode.h"
21 #include "SparseMatrix.h"
22 // #include "NeuroStencil.h"
23 #include "NeuroMesh.h"
24 #include "../utility/Vec.h"
25 #include "CylMesh.h"
26 #include "SpineEntry.h"
27 #include "SpineMesh.h"
28 #include "PsdMesh.h"
29 
34 void testVolScaling()
35 {
36  extern Id makeReacTest();
37  Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
38 
39  Id kin = makeReacTest();
40  Id subCompt = s->doCreate( "CubeMesh", kin, "subCompt", 1 );
41  Field< double >::set( subCompt, "volume", 1e-16 );
42  Id SP = s->doCreate( "Pool", subCompt, "SP", 1 );
43  Field< double >::set( SP, "concInit", 2 );
44 
45  vector< double > n;
46  n.push_back( Field< double >::get( ObjId( "/kinetics/A" ), "nInit" ) );
47  n.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool" ), "nInit" ) );
48  n.push_back( Field< double >::get( ObjId( "/kinetics/r1" ), "numKf" ));
49  n.push_back( Field< double >::get( ObjId( "/kinetics/r1" ), "numKb" ));
50  n.push_back( Field< double >::get( ObjId( "/kinetics/r2" ), "numKf" ));
51  n.push_back( Field< double >::get( ObjId( "/kinetics/r2" ), "numKb" ));
52  n.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k1" ) );
53  n.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k2" ) );
54  n.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k3" ) );
55  n.push_back( Field< double >::get( ObjId( "/kinetics/e2Pool/e2" ), "Km" ) );
56  n.push_back( Field< double >::get( ObjId( "/kinetics/e2Pool/e2" ), "kcat" ) );
57  n.push_back( Field< double >::get( SP, "nInit" ) );
58 
59  double vol = Field< double >::get( kin, "volume" );
60 
61  vol *= 10;
62  Field< double >::set( kin, "volume", vol );
63 
64  vector< double > m;
65  m.push_back( Field< double >::get( ObjId( "/kinetics/A" ), "nInit" ) );
66  m.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool" ), "nInit" ) );
67  m.push_back( Field< double >::get( ObjId( "/kinetics/r1" ), "numKf" ));
68  m.push_back( Field< double >::get( ObjId( "/kinetics/r1" ), "numKb" ));
69  m.push_back( Field< double >::get( ObjId( "/kinetics/r2" ), "numKf" ));
70  m.push_back( Field< double >::get( ObjId( "/kinetics/r2" ), "numKb" ));
71  m.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k1" ) );
72  m.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k2" ) );
73  m.push_back( Field< double >::get( ObjId( "/kinetics/e1Pool/e1" ), "k3" ) );
74  m.push_back( Field< double >::get( ObjId( "/kinetics/e2Pool/e2" ), "Km" ) );
75  m.push_back( Field< double >::get( ObjId( "/kinetics/e2Pool/e2" ), "kcat" ) );
76  m.push_back( Field< double >::get( SP, "nInit" ) );
77 
78  assert( doubleEq( n[0] * 10, m[0] ) ); // A nInit
79  assert( doubleEq( n[1] * 10, m[1] ) ); // e1Pool
80  assert( doubleEq( n[2] / 10, m[2] ) ); // r1 numKf
81  assert( doubleEq( n[3], m[3] ) ); // r1 numKb
82  assert( doubleEq( n[4] / 10, m[4] ) ); // r2 numKf
83  assert( doubleEq( n[5], m[5] ) ); // r2 numKb
84  assert( doubleEq( n[6] / 10, m[6] ) ); // e1 k1
85  assert( doubleEq( n[7], m[7] ) ); // e1 k2
86  assert( doubleEq( n[8], m[8] ) ); // e1 k3
87  assert( doubleEq( n[9], m[9] ) ); // e2 Km
88  assert( doubleEq( n[10], m[10] ) ); // e2 kcat
89  assert( doubleEq( n[11], m[11] ) ); // SP nInit
90 
91  s->doDelete( kin );
92  cout << "." << flush;
93 }
94 
95 #if 0
96 
99 void testCylBase()
100 {
101  // CylBase( x, y, z, dia, len, numDivs );
102  CylBase a( 0, 0, 0, 1, 1, 1 );
103  CylBase b( 1, 2, 3, 1, 2, 10 );
104 
105  assert( doubleEq( b.volume( a ), PI * 0.25 * 2 ) );
106  assert( a.getNumDivs() == 1 );
107  assert( b.getNumDivs() == 10 );
108 
109  for ( unsigned int i = 0; i < 10; ++i ) {
110  assert( doubleEq( b.voxelVolume( a, i ), PI * 0.25 * 2 * 0.1 ) );
111  vector< double > coords = b.getCoordinates( a, i );
112  double x = i;
113  assert( doubleEq( coords[0], x / 10.0 ) );
114  assert( doubleEq( coords[1], x * 2.0 / 10.0 ) );
115  assert( doubleEq( coords[2], x * 3.0 / 10.0 ) );
116  assert( doubleEq( coords[3], (1.0 + x) / 10.0 ) );
117  assert( doubleEq( coords[4], (1.0 + x) * 2.0 / 10.0 ) );
118  assert( doubleEq( coords[5], (1.0 + x) * 3.0 / 10.0 ) );
119  assert( doubleEq( coords[6], 0.5 ) );
120  assert( doubleEq( coords[7], 0.5 ) );
121  assert( doubleEq( coords[8], 0.0 ) );
122  assert( doubleEq( coords[9], 0.0 ) );
123 
124  assert( doubleEq( b.getDiffusionArea( a, i ), PI * 0.25 ) );
125  }
126  b.setDia( 2.0 );
127  assert( doubleEq( b.getDia(), 2.0 ) );
128  assert( doubleEq( b.volume( a ), PI * (2*(0.5*0.5+0.5+1)/3.0) ) );
129  for ( unsigned int i = 0; i < 10; ++i ) {
130  double x = static_cast< double >( i ) / 10.0;
131  double r0 = 0.5 * ( a.getDia() * ( 1.0 - x ) + b.getDia() * x );
132  x = x + 0.1;
133  double r1 = 0.5 * ( a.getDia() * ( 1.0 - x ) + b.getDia() * x );
134  // len = 2, and we have 10 segments of it.
135  double vv = 0.2 * ( r0 * r0 + r0 * r1 + r1 * r1 ) * PI / 3.0;
136  assert( doubleEq( b.voxelVolume( a, i ), vv ) );
137  assert( doubleEq( b.getDiffusionArea( a, i ), PI * r0 * r0 ) );
138  }
139 
140  cout << "." << flush;
141 }
142 
146 void testNeuroNode()
147 {
148  CylBase a( 0, 0, 0, 1, 1, 1 );
149  CylBase b( 1, 2, 3, 2, 2, 10 );
150  CylBase dcb( 1, 2, 3, 2, 2, 0 );
151  vector< unsigned int > twoKids( 2, 2 );
152  twoKids[0] = 2;
153  twoKids[1] = 4;
154  vector< unsigned int > oneKid( 1, 2 );
155  vector< unsigned int > noKids( 0 );
156  NeuroNode na( a, 0, twoKids, 0, Id(), true );
157  NeuroNode ndummy( dcb, 0, oneKid, 1, Id(), false );
158  NeuroNode nb( b, 1, noKids, 1, Id(), false );
159 
160  assert( na.parent() == 0 );
161  assert( na.startFid() == 0 );
162  assert( na.elecCompt() == Id() );
163  assert( na.isDummyNode() == false );
164  assert( na.isSphere() == true );
165  assert( na.isStartNode() == true );
166  assert( na.children().size() == 2 );
167  assert( na.children()[0] == 2 );
168  assert( na.children()[1] == 4 );
169  assert( doubleEq( na.volume( a ), PI * 0.25 ) );
170 
171  assert( ndummy.parent() == 0 );
172  assert( ndummy.startFid() == 1 );
173  assert( ndummy.elecCompt() == Id() );
174  assert( ndummy.isDummyNode() == true );
175  assert( ndummy.isSphere() == false );
176  assert( ndummy.isStartNode() == false );
177  assert( ndummy.children().size() == 1 );
178  assert( ndummy.children()[0] == 2 );
179 
180  assert( nb.parent() == 1 );
181  assert( nb.startFid() == 1 );
182  assert( nb.elecCompt() == Id() );
183  assert( nb.isDummyNode() == false );
184  assert( nb.isSphere() == false );
185  assert( nb.isStartNode() == false );
186  assert( nb.children().size() == 0 );
187  assert( doubleEq( nb.volume( a ), PI * (2*(0.5*0.5+0.5+1)/3.0) ) );
188 
189  cout << "." << flush;
190 }
191 
192 void zebraConcPattern( vector< vector< double > >& S )
193 {
194  // Set up an alternating pattern of zero/nonzero voxels to do test
195  // calculations
196  assert( S.size() == 44 );
197  for ( unsigned int i = 0; i < S.size(); ++i )
198  S[i][0] = 0;
199 
200  S[0][0] = 100; // soma
201  S[2][0] = 10; // basal dend
202  S[4][0] = 1; // basal
203  S[6][0] = 10; // basal
204  S[8][0] = 1; // basal
205  S[10][0] = 10; // basal
206  S[12][0] = 1; // basal
207  S[14][0] = 10; // basal
208  S[16][0] = 1; // basal
209  S[18][0] = 10; // basal
210  S[20][0] = 1; // basal tip
211 
212  S[22][0] = 10; // Apical left dend second compt.
213  S[24][0] = 1; // Apical left dend
214  S[26][0] = 10; // Apical left dend
215  S[28][0] = 1; // Apical left dend
216  S[30][0] = 10; // Apical left dend last compt
217 
218  // Tip compts are 31 and 32:
219  S[32][0] = 1; // Apical middle Tip compartment.
220 
221  S[34][0] = 10; // Apical right dend second compt.
222  S[36][0] = 1; // Apical right dend
223  S[38][0] = 10; // Apical right dend
224  S[40][0] = 1; // Apical right dend
225  S[42][0] = 10; // Apical right dend last compt
226 
227  S[43][0] = 1; // Apical right tip.
228 }
229 
232 void uniformConcPattern( vector< vector< double > >& S,
233  vector< double >& vs )
234 {
235  // Set up uniform concs throughout. This means to scale each # by
236  // voxel size.
237  assert( S.size() == 44 );
238  for ( unsigned int i = 0; i < S.size(); ++i )
239  S[i][0] = 10 * vs[i];
240 }
241 
276 unsigned int buildNode( vector< NeuroNode >& nodes, unsigned int parent,
277  double x, double y, double z, double dia,
278  unsigned int numDivs, bool isDummy, unsigned int startFid )
279 {
280  x *= 1e-6;
281  y *= 1e-6;
282  z *= 1e-6;
283  dia *= 1e-6;
284  CylBase cb( x, y, z, dia, dia/2, numDivs );
285  vector< unsigned int > kids( 0 );
286  if ( nodes.size() == 0 ) { // make soma
287  NeuroNode nn( cb, parent, kids, startFid, Id(), true );
288  nodes.push_back( nn );
289  } else if ( isDummy ) { // make dummy
290  NeuroNode nn( cb, parent, kids, startFid, Id(), false );
291  nodes.push_back( nn );
292  } else {
293  NeuroNode nn( cb, parent, kids, startFid, Id(), false );
294  nodes.push_back( nn );
295  }
296  NeuroNode& paNode = nodes[ parent ];
297  nodes.back().calculateLength( paNode );
298  return startFid + numDivs;
299 }
300 
301 void connectChildNodes( vector< NeuroNode >& nodes )
302 {
303  assert( nodes.size() == 11 );
304  for ( unsigned int i = 1; i < nodes.size(); ++i ) {
305  assert( nodes[i].parent() < nodes.size() );
306  NeuroNode& parent = nodes[ nodes[i].parent() ];
307  parent.addChild( i );
308  }
310  // Check children
311  assert( nodes[0].children().size() == 3 );
312  assert( nodes[0].children()[0] == 1 );
313  assert( nodes[0].children()[1] == 4 );
314  assert( nodes[0].children()[2] == 8 );
315  assert( nodes[1].children().size() == 1 );
316  assert( nodes[1].children()[0] == 2 );
317  assert( nodes[2].children().size() == 1 );
318  assert( nodes[2].children()[0] == 3 );
319  assert( nodes[3].children().size() == 0 );
320  assert( nodes[4].children().size() == 1 );
321  assert( nodes[4].children()[0] == 5 );
322  assert( nodes[5].children().size() == 2 );
323  assert( nodes[5].children()[0] == 6 );
324  assert( nodes[5].children()[1] == 7 );
325  assert( nodes[6].children().size() == 0 );
326  assert( nodes[7].children().size() == 0 );
327  assert( nodes[8].children().size() == 1 );
328  assert( nodes[8].children()[0] == 9 );
329  assert( nodes[9].children().size() == 1 );
330  assert( nodes[9].children()[0] == 10 );
331  assert( nodes[10].children().size() == 0 );
332 }
333 
334 
338 void testCylMesh()
339 {
340  CylMesh cm;
341  assert( cm.getMeshType( 0 ) == CYL );
342  assert( cm.getMeshDimensions( 0 ) == 3 );
343  assert( cm.getDimensions() == 3 );
344 
345  vector< double > coords( 9 );
346  coords[0] = 1; // X0
347  coords[1] = 2; // Y0
348  coords[2] = 3; // Z0
349 
350  coords[3] = 3; // X1
351  coords[4] = 5; // Y1
352  coords[5] = 7; // Z1
353 
354  coords[6] = 1; // R0
355  coords[7] = 2; // R1
356 
357  coords[8] = 1; // DiffLength
358 
359  cm.innerSetCoords( coords );
360 
361  assert( doubleEq( cm.getX0(), 1 ) );
362  assert( doubleEq( cm.getY0(), 2 ) );
363  assert( doubleEq( cm.getZ0(), 3 ) );
364  assert( doubleEq( cm.getR0(), 1 ) );
365 
366  assert( doubleEq( cm.getX1(), 3 ) );
367  assert( doubleEq( cm.getY1(), 5 ) );
368  assert( doubleEq( cm.getZ1(), 7 ) );
369  assert( doubleEq( cm.getR1(), 2 ) );
370 
371  assert( doubleEq( cm.getDiffLength(), sqrt( 29.0 ) / 5.0 ) );
372 
373  cm.setX0( 2 );
374  cm.setY0( 3 );
375  cm.setZ0( 4 );
376  cm.setR0( 2 );
377 
378  cm.setX1( 4 );
379  cm.setY1( 6 );
380  cm.setZ1( 8 );
381  cm.setR1( 3 );
382 
383  cm.setDiffLength( 2.0 );
384 
385  vector< double > temp = cm.getCoords( Id().eref(), 0 );
386  assert( temp.size() == 9 );
387  // Can't check on the last coord as it is DiffLength, it changes.
388  for ( unsigned int i = 0; i < temp.size() - 1; ++i )
389  assert( doubleEq( temp[i], coords[i] + 1 ) );
390 
391  double totLen = sqrt( 29.0 );
392  assert( doubleEq( cm.getTotLength() , totLen ) );
393 
394  cm.setDiffLength( 1.0 );
395  assert( cm.getNumEntries() == 5 );
396  assert( doubleEq( cm.getDiffLength(), totLen / 5 ) );
397 
399  assert( doubleEq( cm.getMeshEntryVolume( 2 ), 2.5 * 2.5 * PI * totLen / 5 ) );
400 
402  // LenSlope/totLen = 0.016 =
403  // 1/numEntries * (r1-r0)/numEntries * 2/(r0+r1) = 1/25 * 1 * 2/5
404  // Here are the fractional positions
405  // part0 = 1/5 - 0.032: end= 0.2 - 0.032
406  // part1 = 1/5 - 0.016: end = 0.4 - 0.048
407  // part2 = 1/5 : end = 0.6 - 0.048
408  // part3 = 1/5 + 0.016 : end = 0.8 - 0.032
409  // part4 = 1/5 + 0.032 : end = 1
410 
411  coords = cm.getCoordinates( 2 );
412  assert( coords.size() == 10 );
413  assert( doubleEq( coords[0], 2 + (0.4 - 0.048) * 2 ) );
414  assert( doubleEq( coords[1], 3 + (0.4 - 0.048) * 3 ) );
415  assert( doubleEq( coords[2], 4 + (0.4 - 0.048) * 4 ) );
416 
417  assert( doubleEq( coords[3], 2 + (0.6 - 0.048) * 2 ) );
418  assert( doubleEq( coords[4], 3 + (0.6 - 0.048) * 3 ) );
419  assert( doubleEq( coords[5], 4 + (0.6 - 0.048) * 4 ) );
420 
421  assert( doubleEq( coords[6], 2.4 ) );
422  assert( doubleEq( coords[7], 2.6 ) );
423 
425  vector< unsigned int > neighbors = cm.getNeighbors( 2 );
426  assert( neighbors.size() == 2 );
427  assert( neighbors[0] == 1 );
428  assert( neighbors[1] == 3 );
429 
431  coords = cm.getDiffusionArea( 2 );
432  assert( coords.size() == 2 );
433  assert( doubleEq( coords[0], 2.4 * 2.4 * PI ) );
434  assert( doubleEq( coords[1], 2.6 * 2.6 * PI ) );
435 
437  coords = cm.getCoords( Id().eref(), 0 );
438  assert( coords.size() == 9 );
439 
440  double x, y, z;
441  double ne = cm.getNumEntries();
442  double ux = ( coords[3] - coords[0] ) / ne;
443  double uy = ( coords[4] - coords[1] ) / ne;
444  double uz = ( coords[5] - coords[2] ) / ne;
445  cm.indexToSpace( 0, x, y, z );
446  assert( doubleEq( x, coords[0] + 0.5 * ux ) );
447  assert( doubleEq( y, coords[1] + 0.5 * uy ) );
448  assert( doubleEq( z, coords[2] + 0.5 * uz ) );
449  cm.indexToSpace( 4, x, y, z );
450  assert( doubleEq( x, coords[3] - 0.5 * ux ) );
451  assert( doubleEq( y, coords[4] - 0.5 * uy ) );
452  assert( doubleEq( z, coords[5] - 0.5 * uz ) );
453 
455 
456  unsigned int index;
457  double dist = cm.nearest( x, y, z, index );
458  assert( doubleEq( dist, 0.0 ) );
459  assert( index == cm.innerGetNumEntries() - 1 );
460  x += 1; // Should be within R1
461  dist = cm.nearest( x, y, z, index );
462  assert( dist > 0.0 && dist < 1.0 );
463  assert( index == cm.innerGetNumEntries() - 1 );
464 
465  x += 5; // Should now be outside R1
466  dist = cm.nearest( x, y, z, index );
467  assert( dist < 0.0 );
468  assert( index == cm.innerGetNumEntries() - 1 );
469 
471  dist = cm.selectGridVolume( 10.0 );
472  assert( doubleEq( dist, 0.1 * cm.getDiffLength() ) );
473  dist = cm.selectGridVolume( 0.1 );
474  assert( dist <= 0.01 );
475 
477  // Test CylMesh::fillPointsOnCircle
478  CubeMesh cube;
479  cube.setPreserveNumEntries( 0 );
480  coords[0] = 0; // X0
481  coords[1] = 0; // Y0
482  coords[2] = 0; // Z0
483 
484  coords[3] = 100; // X1
485  coords[4] = 100; // Y1
486  coords[5] = 100; // Z1
487  coords[6] = 100; // DX
488  coords[7] = 100; // DY
489  coords[8] = 100; // DZ
490  cube.innerSetCoords( coords );
491  vector< VoxelJunction > ret;
492  cm.matchCubeMeshEntries( &cube, ret );
493  assert( ret.size() == cm.innerGetNumEntries() );
494  for ( unsigned int i = 0; i < ret.size(); ++i ) {
495  assert( ret[i].first == i );
496  assert( ret[i].second == 0 ); // Only one cube entry
497  double r = cm.getR0() +
498  ( cm.getR1() - cm.getR0() ) *
499  cm.getDiffLength() * ( 0.5 + i ) / cm.getTotLength();
500  double a = r * cm.getDiffLength() * 2 * PI;
501  //assert( doubleApprox( ret[i].diffScale, a ) );
502  // cout << i << ". mesh: " << ret[i].diffScale << ", calc: " << a << endl;
503  assert( fabs( ret[i].diffScale - a ) < 0.5 );
504  }
505 
507  // We're going to set up a new cylinder to abut the old one. The
508  // coords of x1 of the new cylinder are the same x0 of the old one, but
509  // the coords of the other end change.
510  coords = cm.getCoords( Id().eref(), 0 );
511  coords[3] = coords[0];
512  coords[4] = coords[1];
513  coords[5] = coords[2];
514  coords[2] -= 10; // x and y stay put, it is a vertical cylinder.
515  coords[6] = 1; // R0
516  coords[7] = 1; // R1
517  coords[8] = 1; // DiffLength
518  CylMesh cm2;
519  cm2.innerSetCoords( coords );
520  assert( doubleEq( cm2.getDiffLength(), 1 ) );
521  assert( cm2.getNumEntries() == 10 );
522  vector< VoxelJunction > vj;
523  cm.matchMeshEntries( &cm2, vj );
524  assert( vj.size() == 1 );
525  assert( vj[0].first == 0 );
526  assert( vj[0].second == cm2.getNumEntries() - 1 );
527  double xda = 2.0 * 1 * 1 * PI / ( cm.getDiffLength() + cm2.getDiffLength() );
528  assert( doubleEq( vj[0].diffScale, xda ) );
529 
530  cout << "." << flush;
531 }
532 
533 
537 void testMidLevelCylMesh()
538 {
539  Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
540 
541  Id cylId = s->doCreate( "CylMesh", Id(), "cyl", 1 );
542  Id meshId( cylId.value() + 1 );
543 
544  vector< double > coords( 9 );
545  coords[0] = 1; // X0
546  coords[1] = 2; // Y0
547  coords[2] = 3; // Z0
548 
549  coords[3] = 3; // X1
550  coords[4] = 5; // Y1
551  coords[5] = 7; // Z1
552 
553  coords[6] = 1; // R0
554  coords[7] = 2; // R1
555 
556  coords[8] = 1; // DiffLength
557 
558  bool ret = Field< vector< double > >::set( cylId, "coords", coords );
559  assert( ret );
560 
561  assert( doubleEq( Field< double >::get( cylId, "x0" ), 1 ) );
562  assert( doubleEq( Field< double >::get( cylId, "y0" ), 2 ) );
563  assert( doubleEq( Field< double >::get( cylId, "z0" ), 3 ) );
564  assert( doubleEq( Field< double >::get( cylId, "x1" ), 3 ) );
565  assert( doubleEq( Field< double >::get( cylId, "y1" ), 5 ) );
566  assert( doubleEq( Field< double >::get( cylId, "z1" ), 7 ) );
567  assert( doubleEq( Field< double >::get( cylId, "r0" ), 1 ) );
568  assert( doubleEq( Field< double >::get( cylId, "r1" ), 2 ) );
569 
570  ret = Field< double >::set( cylId, "diffLength", 1 );
571  assert( ret );
572  meshId.element()->syncFieldDim();
573 
574  assert( meshId()->dataHandler()->localEntries() == 5 );
575 
576  unsigned int n = Field< unsigned int >::get( cylId, "num_mesh" );
577  assert( n == 5 );
578 
579  ObjId oid( meshId, DataId( 2 ) );
580 
581  double totLen = sqrt( 29.0 );
582  assert( doubleEq( Field< double >::get( oid, "volume" ),
583  1.5 * 1.5 * PI * totLen / 5 ) );
584 
585  vector< unsigned int > neighbors =
586  Field< vector< unsigned int > >::get( oid, "neighbors" );
587  assert( neighbors.size() == 2 );
588  assert( neighbors[0] = 1 );
589  assert( neighbors[1] = 3 );
590 
591  s->doDelete( cylId );
592 
593  cout << "." << flush;
594 }
595 #endif
596 
600 void testCubeMesh()
601 {
602  CubeMesh cm;
603  cm.setPreserveNumEntries( 0 );
604  assert( cm.getMeshType( 0 ) == CUBOID );
605  assert( cm.getMeshDimensions( 0 ) == 3 );
606  assert( cm.getDimensions() == 3 );
607 
608  vector< double > coords( 9 );
609  coords[0] = 0; // X0
610  coords[1] = 0; // Y0
611  coords[2] = 0; // Z0
612 
613  coords[3] = 2; // X1
614  coords[4] = 4; // Y1
615  coords[5] = 8; // Z1
616 
617  coords[6] = 1; // DX
618  coords[7] = 1; // DY
619  coords[8] = 1; // DZ
620 
621  cm.innerSetCoords( coords );
622 
623  vector< unsigned int > neighbors = cm.getNeighbors( 0 );
624  assert( neighbors.size() == 3 );
625  assert( neighbors[0] = 1 );
626  assert( neighbors[0] = 2 );
627  assert( neighbors[0] = 8 );
628 
629  assert( cm.innerGetNumEntries() == 64 );
630  assert( doubleEq( cm.getX0(), 0 ) );
631  assert( doubleEq( cm.getY0(), 0 ) );
632  assert( doubleEq( cm.getZ0(), 0 ) );
633 
634  assert( doubleEq( cm.getX1(), 2 ) );
635  assert( doubleEq( cm.getY1(), 4 ) );
636  assert( doubleEq( cm.getZ1(), 8 ) );
637 
638  assert( doubleEq( cm.getDx(), 1 ) );
639  assert( doubleEq( cm.getDy(), 1 ) );
640  assert( doubleEq( cm.getDz(), 1 ) );
641 
642  assert( cm.getNx() == 2 );
643  assert( cm.getNy() == 4 );
644  assert( cm.getNz() == 8 );
645 
646  cm.setX0( 1 );
647  cm.setY0( 2 );
648  cm.setZ0( 4 );
649 
650  cm.setX1( 5 );
651  cm.setY1( 6 );
652  cm.setZ1( 8 );
653 
654  vector< double > temp = cm.getCoords( Id().eref() );
655  assert( temp.size() == 9 );
656  assert( doubleEq( temp[0], 1 ) );
657  assert( doubleEq( temp[1], 2 ) );
658  assert( doubleEq( temp[2], 4 ) );
659  assert( doubleEq( temp[3], 5 ) );
660  assert( doubleEq( temp[4], 6 ) );
661  assert( doubleEq( temp[5], 8 ) );
662  assert( doubleEq( temp[6], 1 ) );
663  assert( doubleEq( temp[7], 1 ) );
664  assert( doubleEq( temp[8], 1 ) );
665  assert( cm.innerGetNumEntries() == 64 );
666  assert( cm.getNx() == 4 );
667  assert( cm.getNy() == 4 );
668  assert( cm.getNz() == 4 );
669 
670  neighbors = cm.getNeighbors( 0 );
671  assert( neighbors.size() == 3 );
672  assert( neighbors[0] == 1 );
673  assert( neighbors[1] == 4 );
674  assert( neighbors[2] = 16 );
675 
676  neighbors = cm.getNeighbors( 63 );
677  assert( neighbors.size() == 3 );
678  assert( neighbors[0] == 47 );
679  assert( neighbors[1] == 59 );
680  assert( neighbors[2] == 62 );
681 
682  neighbors = cm.getNeighbors( 2 );
683  assert( neighbors.size() == 4 );
684  assert( neighbors[0] == 1 );
685  assert( neighbors[1] == 3 );
686  assert( neighbors[2] == 6 );
687  assert( neighbors[3] == 18 );
688 
689  neighbors = cm.getNeighbors( 6 );
690  assert( neighbors.size() == 5 );
691  assert( neighbors[0] == 2 );
692  assert( neighbors[1] == 5 );
693  assert( neighbors[2] == 7 );
694  assert( neighbors[3] == 10 );
695  assert( neighbors[4] == 22 );
696 
697  neighbors = cm.getNeighbors( 22 );
698  assert( neighbors.size() == 6 );
699  assert( neighbors[0] == 6 );
700  assert( neighbors[1] == 18 );
701  assert( neighbors[2] == 21 );
702  assert( neighbors[3] == 23 );
703  assert( neighbors[4] == 26 );
704  assert( neighbors[5] == 38 );
705 
706  cm.setPreserveNumEntries( 1 );
707  assert( cm.getNx() == 4 );
708  assert( cm.getNy() == 4 );
709  assert( cm.getNz() == 4 );
710  assert( doubleEq( cm.getDx(), 1.0 ) );
711  assert( doubleEq( cm.getDy(), 1.0 ) );
712  assert( doubleEq( cm.getDz(), 1.0 ) );
713 
714  cm.setX0( 0 );
715  cm.setY0( 0 );
716  cm.setZ0( 0 );
717  // x1 is 5, y1 is 6 and z1 is 8
718 
719  assert( doubleEq( cm.getDx(), 1.25 ) );
720  assert( doubleEq( cm.getDy(), 1.5 ) );
721  assert( doubleEq( cm.getDz(), 2.0 ) );
722 
723  cout << "." << flush;
724 }
725 
726 void testCubeMeshExtendStencil()
727 {
728  CubeMesh cm0;
729  cm0.setPreserveNumEntries( 0 );
730  assert( cm0.getMeshType( 0 ) == CUBOID );
731  assert( cm0.getMeshDimensions( 0 ) == 3 );
732  assert( cm0.getDimensions() == 3 );
733  CubeMesh cm1 = cm0;
734 
735  vector< double > coords( 9 );
736  coords[0] = 0; // X0
737  coords[1] = 0; // Y0
738  coords[2] = 0; // Z0
739 
740  coords[3] = 2; // X1
741  coords[4] = 4; // Y1
742  coords[5] = 8; // Z1
743 
744  coords[6] = 1; // DX
745  coords[7] = 1; // DY
746  coords[8] = 1; // DZ
747 
748  cm0.innerSetCoords( coords );
749  coords[2] = 8; // 2x4 face abuts.
750  coords[5] = 16;
751  cm1.innerSetCoords( coords );
752 
753  const double* entry;
754  const unsigned int* colIndex;
755  unsigned int num = cm0.getStencilRow( 0, &entry, &colIndex );
756  assert( num == 3 );
757  assert( colIndex[0] == 1 );
758  assert( colIndex[1] == 2 );
759  assert( colIndex[2] == 8 );
760 
761  num = cm0.getStencilRow( 56, &entry, &colIndex );
762  assert( num == 3 );
763  assert( colIndex[0] == 48 );
764  assert( colIndex[1] == 57 );
765  assert( colIndex[2] == 58 );
766 
767  vector< VoxelJunction > vj;
768  for ( unsigned int i = 0; i < 8; ++i ) {
769  vj.push_back( VoxelJunction( 56 + i, i ) );
770  }
771  cm0.extendStencil( &cm1, vj );
772 
773  num = cm0.getStencilRow( 56, &entry, &colIndex );
774  assert( num == 4 );
775  assert( colIndex[0] == 48 );
776  assert( colIndex[1] == 57 );
777  assert( colIndex[2] == 58 );
778  assert( colIndex[3] == 64 );
779 
780  for ( unsigned int i = 0; i < 8; ++i ) {
781  num = cm0.getStencilRow( 64 + i, &entry, &colIndex );
782  assert( num == 1 );
783  assert( colIndex[0] == 56 + i );
784  }
785 
786  cout << "." << flush;
787 }
788 
792 void testReMesh()
793 {
794  Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
795  Id base = s->doCreate( "Neutral", Id(), "base", 1 );
796 
797  Id cube = s->doCreate( "CubeMesh", base, "cube", 1 );
798  bool ret = SetGet2< double, unsigned int >::set(
799  cube, "buildDefaultMesh", 1.0, 1 );
800  assert( ret );
801  unsigned int vol = Field< double >::get( cube, "volume" );
802  assert( doubleEq( vol, 1.0 ) );
803  Id pool = s->doCreate( "Pool", cube, "pool", 1 );
804  Id mesh( "/base/cube/mesh" );
805  assert( mesh != Id() );
806 
808  // 1 millimolar in 1 m^3 is 1 mole per liter.
809  double linsize = Field< double >::get( pool, "volume" );
810  assert( doubleEq( linsize, 1.0 ) );
811 
812  ret = Field< double >::set( pool, "conc", 1 );
813  assert( ret );
814  double n = Field< double >::get( pool, "n" );
815  assert( doubleEq( n, NA ) );
816 
817  ret = SetGet2< double, unsigned int >::set(
818  cube, "buildDefaultMesh", 1.0e-3, 1 );
819  Field< double >::set( pool, "conc", 1 );
820  n = Field< double >::get( pool, "n" );
821  assert( doubleEq( n, NA / 1000.0 ) );
822 
824  // Next we do the remeshing.
825  double x = 1.234;
826  Field< double >::set( pool, "concInit", x );
827  ret = SetGet2< double, unsigned int >::set(
828  cube, "buildDefaultMesh", 1, 8 );
829  // This is nasty, needs the change to propagate through messages.
830  linsize = Field< double >::get( pool, "volume" );
831  assert( doubleEq( linsize, 0.125 ) );
832 
833  n = Field< double >::get( ObjId( pool, 0 ), "concInit" );
834  assert( doubleEq( n, x ) );
835  n = Field< double >::get( ObjId( pool, 7 ), "concInit" );
836  assert( doubleEq( n, x ) );
837  n = Field< double >::get( ObjId( pool, 0 ), "nInit" );
838  assert( doubleEq( n, x * NA / 8.0 ) );
839  n = Field< double >::get( ObjId( pool, 7 ), "nInit" );
840  assert( doubleEq( n, x * NA / 8.0 ) );
841  n = Field< double >::get( ObjId( pool, 0 ), "conc" );
842  assert( doubleEq( n, x ) );
843  n = Field< double >::get( ObjId( pool, 7 ), "conc" );
844  assert( doubleEq( n, x ) );
845 
847  s->doDelete( base );
848  cout << "." << flush;
849 }
850 
855 Id makeCompt( Id parentCompt, Id parentObj,
856  string name, double len, double dia, double theta )
857 {
858  Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
859  Id ret = shell->doCreate( "Compartment", parentObj, name, 1 );
860  double pax = 0;
861  double pay = 0;
862  if ( parentCompt != Id() ) {
863  pax = Field< double >::get( parentCompt, "x" );
864  pay = Field< double >::get( parentCompt, "y" );
865  shell->doAddMsg( "Single", parentCompt, "raxial", ret, "axial" );
866  }
867  Field< double >::set( ret, "x0", pax );
868  Field< double >::set( ret, "y0", pay );
869  Field< double >::set( ret, "z0", 0.0 );
870  double x = pax + len * cos( theta * PI / 180.0 );
871  double y = pay + len * sin( theta * PI / 180.0 );
872  Field< double >::set( ret, "x", x );
873  Field< double >::set( ret, "y", y );
874  Field< double >::set( ret, "z", 0.0 );
875  Field< double >::set( ret, "diameter", dia );
876  Field< double >::set( ret, "length", len );
877 
878  return ret;
879 }
880 
881 #if 0
882 // Assumes parent dend is along x axis.
883 Id makeSpine( Id parentCompt, Id parentObj, unsigned int index,
884  double frac, double len, double dia, double theta )
885 {
886  const double RA = 1;
887  const double RM = 1;
888  const double CM = 0.01;
889  assert( parentCompt != Id() );
890  Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
891  double pax0 = Field< double >::get( parentCompt, "x0" );
892  double pay0 = Field< double >::get( parentCompt, "y0" );
893  double paz0 = Field< double >::get( parentCompt, "z0" );
894  double pax1 = Field< double >::get( parentCompt, "x" );
895  double pay1 = Field< double >::get( parentCompt, "y" );
896  double paz1 = Field< double >::get( parentCompt, "z" );
897 
898  stringstream ss;
899  ss << "shaft" << index;
900  string sname = ss.str();
901  stringstream ss2;
902  ss2 << "head" << index;
903  string hname = ss2.str();
904 
905  Id shaft = shell->doCreate( "Compartment", parentObj, sname, 1 );
906  ObjId mid =
907  shell->doAddMsg( "Single", parentCompt, "raxial", shaft, "axial" );
908  assert( !mid.bad() );
909  double x = pax0 + frac * ( pax1 - pax0 );
910  double y = pay0 + frac * ( pay1 - pay0 );
911  double z = paz0 + frac * ( paz1 - paz0 );
912  Field< double >::set( shaft, "x0", x );
913  Field< double >::set( shaft, "y0", y );
914  Field< double >::set( shaft, "z0", z );
915  double sy = y + len * cos( theta * PI / 180.0 );
916  double sz = z + len * sin( theta * PI / 180.0 );
917  Field< double >::set( shaft, "x", x );
918  Field< double >::set( shaft, "y", sy );
919  Field< double >::set( shaft, "z", sz );
920  Field< double >::set( shaft, "diameter", dia/10.0 );
921  Field< double >::set( shaft, "length", len );
922  double xa = PI * dia * dia / 400.0;
923  double circumference = PI * dia / 10.0;
924  Field< double >::set( shaft, "Ra", RA * len / xa );
925  Field< double >::set( shaft, "Rm", RM / ( len * circumference ) );
926  Field< double >::set( shaft, "Cm", CM * len * circumference );
927 
928  Id head = shell->doCreate( "Compartment", parentObj, hname, 1 );
929  mid = shell->doAddMsg( "Single", shaft, "raxial", head, "axial" );
930  assert( !mid.bad() );
931  Field< double >::set( head, "x0", x );
932  Field< double >::set( head, "y0", sy );
933  Field< double >::set( head, "z0", sz );
934  double hy = sy + len * cos( theta * PI / 180.0 );
935  double hz = sz + len * sin( theta * PI / 180.0 );
936  Field< double >::set( head, "x", x );
937  Field< double >::set( head, "y", hy );
938  Field< double >::set( head, "z", hz );
939  Field< double >::set( head, "diameter", dia );
940  Field< double >::set( head, "length", len );
941  xa = PI * dia * dia / 4.0;
942  circumference = PI * dia;
943  Field< double >::set( head, "Ra", RA * len / xa );
944  Field< double >::set( head, "Rm", RM / ( len * circumference ) );
945  Field< double >::set( head, "Cm", CM * len * circumference );
946 
947  return head;
948 }
949 
964 pair< unsigned int, unsigned int > buildBranchingCell(
965  Id cell, double len, double dia )
966 {
967  Id soma = makeCompt( Id(), cell, "soma", dia, dia, 90 );
968  dia /= sqrt( 2.0 );
969  Id d1 = makeCompt( soma, cell, "d1", len , dia, 0 );
970  Id d2 = makeCompt( soma, cell, "d2", len , dia, 180 );
971  Id d1a = makeCompt( d1, cell, "d1a", len , dia, 0 );
972  Id d2a = makeCompt( d2, cell, "d2a", len , dia, 180 );
973  dia /= sqrt( 2.0 );
974  Id d11 = makeCompt( d1a, cell, "d11", len , dia, -45 );
975  Id d12 = makeCompt( d1a, cell, "d12", len , dia, 45 );
976  Id d21 = makeCompt( d2a, cell, "d21", len , dia, 45 );
977  Id d22 = makeCompt( d2a, cell, "d22", len , dia, -45 );
978  dia /= sqrt( 2.0 );
979  Id d111 = makeCompt( d11, cell, "d111", len , dia, -90 );
980  Id d112 = makeCompt( d11, cell, "d112", len , dia, 0 );
981  Id d121 = makeCompt( d12, cell, "d121", len , dia, 0 );
982  Id d122 = makeCompt( d12, cell, "d122", len , dia, 90 );
983  Id d211 = makeCompt( d21, cell, "d211", len , dia, 90 );
984  Id d212 = makeCompt( d21, cell, "d212", len , dia, 180 );
985  Id d221 = makeCompt( d22, cell, "d221", len , dia, 180 );
986  Id d222 = makeCompt( d22, cell, "d222", len , dia, -90 );
987 
988  return pair< unsigned int, unsigned int >( 17, 161 );
989 }
990 
991  // y = initConc * dx * (0.5 / sqrt( PI * DiffConst * runtime ) ) *
992  // exp( -x * x / ( 4 * DiffConst * runtime ) )
993 double diffusionFunction( double D, double dx, double x, double t )
994 {
995  return
996  dx * (0.5 / sqrt( PI * D * t ) ) * exp( -x * x / ( 4 * D * t ) );
997 }
998 
999 void testNeuroMeshLinear()
1000 {
1001  Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
1002  // Build a linear cylindrical cell, no tapering.
1003  Id cell = shell->doCreate( "Neutral", Id(), "cell", 1 );
1004  unsigned int numCompts = 500;
1005  double dia = 1e-6; // metres
1006  double diffLength = 0.2e-6; // metres
1007  double len = diffLength * numCompts;
1008  double D = 1e-12;
1009  double totNum = 1e6;
1010 
1011  Id soma = makeCompt( Id(), cell, "soma", len, dia, 0 );
1012 
1013  // Scan it with neuroMesh and check outcome.
1014  Id nm = shell->doCreate( "NeuroMesh", Id(), "neuromesh", 1 );
1015  Field< double >::set( nm, "diffLength", diffLength );
1016  Field< string >::set( nm, "geometryPolicy", "cylinder" );
1017  Field< Id >::set( nm, "cell", cell );
1018  unsigned int ns = Field< unsigned int >::get( nm, "numSegments" );
1019  assert( ns == 1 );
1020  unsigned int ndc = Field< unsigned int >::get( nm, "numDiffCompts" );
1021  assert( ndc == numCompts );
1022  const vector< NeuroNode >& nodes =
1023  reinterpret_cast< NeuroMesh* >( nm.eref().data() )->
1024  getNodes();
1025  assert( nodes.size() == 2 ); // Self plus dummy parent.
1026  assert( nodes[0].children().size() == 0 );
1027 
1028  // Insert a molecule at first subdivision of soma. I use a dummy
1029  // matrix S rather than the one in the system.
1030  // Field< double >::set( ObjId( soma, 0 ), "nInit", 1.0e6 );
1031  vector< double > molNum( 1, 0 ); // We only have one pool
1032  // S[meshIndex][poolIndex]
1033  vector< vector< double > > S( ndc, molNum );
1034  S[0][0] = totNum;
1035  vector< double > diffConst( 1, D );
1036  vector< double > temp( 1, 0.0 );
1037  vector< vector< double > > flux( ndc, temp );
1038 
1039  MeshCompt* mc = reinterpret_cast< MeshCompt* >( nm.eref().data() );
1040  assert( mc->getNumEntries() == numCompts );
1041  const double adx = dia * dia * PI * 0.25 / diffLength;
1042  for ( unsigned int i = 0; i < numCompts; ++i ) {
1043  const double* entry;
1044  const unsigned int* colIndex;
1045  unsigned int numAbut = mc->getStencilRow( i, &entry, &colIndex );
1046  if ( i == 0 ) {
1047  assert( numAbut == 1 );
1048  assert( doubleEq( entry[0], adx ) );
1049  assert( colIndex[0] == 1 );
1050  } else if ( i == numCompts - 1 ) {
1051  assert( numAbut == 1 );
1052  assert( doubleEq( entry[0], adx ) );
1053  assert( colIndex[0] == numCompts - 2 );
1054  } else {
1055  assert( numAbut == 2 );
1056  assert( doubleEq( entry[0], adx ) );
1057  assert( colIndex[0] == i - 1 );
1058  assert( doubleEq( entry[1], adx ) );
1059  assert( colIndex[1] == i + 1 );
1060  }
1061  }
1062 
1064  // Test the 'nearest' function
1065  unsigned int index;
1066  double near = mc->nearest( diffLength * 27.5, diffLength / 10.0, 0, index );
1067  assert( index == 27 );
1068  assert( doubleEq( near, diffLength / 10.0 ) );
1069 
1070  near = mc->nearest( -10, 0, 0, index );
1071  assert( index == 0 );
1072  assert( doubleEq( near, -1.0 ) );
1073 
1074  // Test indexToSpace
1075  double x, y, z;
1076  mc->indexToSpace( 27, x, y, z );
1077  assert( doubleEq( x, diffLength * 27.5 ) );
1078  assert( doubleEq( y, 0.0 ) );
1079  assert( doubleEq( z, 0.0 ) );
1080 
1082  shell->doDelete( cell );
1083  shell->doDelete( nm );
1084  cout << "." << flush;
1085 }
1086 
1087 void testNeuroMeshBranching()
1088 {
1089  Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
1090  // Build a cell
1091  Id cell = shell->doCreate( "Neutral", Id(), "cell", 1 );
1092  double len = 10e-6; // metres
1093  double dia = 1e-6; // metres
1094  double diffLength = 1e-6; // metres
1095  double D = 4e-12; // Diff const, m^2/s
1096  double totNum = 1e6; // Molecules
1097  double maxt = 10.0;
1098  double dt = 0.001;
1099 
1100  pair< unsigned int, unsigned int > ret =
1101  buildBranchingCell( cell, len, dia );
1102 
1103  // Scan it with neuroMesh and check outcome.
1104  Id nm = shell->doCreate( "NeuroMesh", Id(), "neuromesh", 1 );
1105  Field< double >::set( nm, "diffLength", diffLength );
1106  Field< string >::set( nm, "geometryPolicy", "cylinder" );
1107  Field< Id >::set( nm, "cell", cell );
1108  unsigned int ns = Field< unsigned int >::get( nm, "numSegments" );
1109  assert( ns == ret.first );
1110  unsigned int ndc = Field< unsigned int >::get( nm, "numDiffCompts" );
1111  assert( ndc == ret.second );
1112  NeuroMesh* neuro = reinterpret_cast< NeuroMesh* >( nm.eref().data() );
1113  const vector< NeuroNode >& nodes = neuro-> getNodes();
1114  assert( nodes.size() == ns + 15 ); // 15 dummy nodes.
1115  // We have depth-first ordering in the nodes vector.
1130  assert( nodes[0].children().size() == 2 ); // soma
1131  assert( nodes[1].children().size() == 1 ); // d1
1132  assert( nodes[2].children().size() == 2 ); // d1a
1133  assert( nodes[3].children().size() == 2 ); // d11
1134  assert( nodes[4].children().size() == 0 ); // d111
1135  assert( nodes[5].children().size() == 0 ); // d112
1136  assert( nodes[6].children().size() == 2 ); // d12
1137  assert( nodes[7].children().size() == 0 ); // d121
1138  assert( nodes[8].children().size() == 0 ); // d122
1139 
1140  assert( nodes[9].children().size() == 1 ); // d2
1141  assert( nodes[10].children().size() == 2 ); // d2a
1142  assert( nodes[11].children().size() == 2 ); // d21
1143  assert( nodes[12].children().size() == 0 ); // d211
1144  assert( nodes[13].children().size() == 0 ); // d212
1145  assert( nodes[14].children().size() == 2 ); // d22
1146  assert( nodes[15].children().size() == 0 ); // d221
1147  assert( nodes[16].children().size() == 0 ); // d222
1148 
1149  // Insert a molecule at soma
1150  vector< double > molNum( 1, 0 ); // We only have one pool
1151  // S[meshIndex][poolIndex]
1152  vector< double > S( ndc, 0.0 );
1153  S[0] = totNum;
1154  vector< double > flux( ndc, 0.0 );
1155  vector< double > vol( ndc, 0.0 );
1156  for ( unsigned int i = 0; i < ndc; ++i )
1157  vol[i] = neuro->getMeshEntryVolume( i );
1158 
1159  assert( doubleEq( vol[0], dia * dia * 0.25 * PI * diffLength ) );
1160  assert( doubleEq( vol[1], dia * dia * 0.125 * PI * diffLength ) );
1161  // Watch diffusion using stencil and direct calls to the flux
1162  // calculations rather than going through the ksolve.
1163  for ( double t = 0; t < maxt; t += dt ) {
1164  flux.assign( ndc, 0.0 );
1165 
1166  for ( unsigned int i = 0; i < ndc; ++i ) {
1167  const double* entry;
1168  const unsigned int* colIndex;
1169  unsigned int numEntries =
1170  neuro->getStencilRow( i, &entry, &colIndex);
1171  for ( unsigned int j = 0; j < numEntries; ++j ) {
1172  unsigned int k = colIndex[j];
1173  double delta = ( S[k]/vol[k] - S[i]/vol[i] ) * entry[j];
1174  flux[ k ] -= delta;
1175  // flux[ i ] += delta; // Am I double counting?
1176  }
1177  }
1178  for ( unsigned int i = 0; i < ndc; ++i )
1179  S[i] += flux[i] * D * dt;
1180  }
1181  double tot = 0.0;
1182  for ( unsigned int i = 0; i < ndc; ++i ) {
1183  tot += S[i];
1184  }
1185 
1186  // Compare with analytic solution.
1187  assert( doubleEq( tot, totNum ) );
1188  double x = 1.5 * diffLength;
1189  for ( unsigned int i = 1; i < 11; ++i ) { // First segment of tree.
1190  double y = totNum * diffusionFunction( D, diffLength, x, maxt );
1191  // cout << "S[" << i << "] = " << S[i] << " " << y << endl;
1192  assert( doubleApprox( y, S[i] ) );
1193  x += diffLength;
1194  }
1195  tot = 0;
1196  for ( double x = diffLength / 2.0 ; x < 10 * len; x += diffLength ) {
1197  // the 2x is needed because we add up only half of the 2-sided
1198  // diffusion distribution.
1199  double y = 2 * totNum * diffusionFunction( D, diffLength, x, maxt);
1200  // cout << floor( x / diffLength ) << ": " << y << endl;
1201  tot += y;
1202  }
1203  assert( doubleEq( tot, totNum ) );
1204 
1206  // Test the 'nearest' function
1207  MeshCompt* mc = reinterpret_cast< MeshCompt* >( nm.eref().data() );
1208  assert( mc->getNumEntries() == ndc );
1209  unsigned int index;
1210  double y, z;
1211  mc->indexToSpace( 60, x, y, z );
1212 
1213  double near = mc->nearest( x, y, z + diffLength/10.0, index );
1214  assert( index == 60 );
1215  assert( doubleEq( near, diffLength / 10.0 ) );
1216 
1217  near = mc->nearest( -10, 0, 0, index );
1218  assert( index == 0 );
1219  assert( doubleEq( near, -1.0 ) );
1220 
1221  vector< NeuroNode > nn = dynamic_cast< NeuroMesh* >( mc )->getNodes();
1222  assert( nn.size() == 32 ); // 17 plus 15 dummies
1223  /*
1224  for ( unsigned int i = 0; i < nn.size(); ++i ) {
1225  if ( !nn[i].isDummyNode() ) {
1226  cout << i << " (" <<
1227  nn[i].getX() << ", " << nn[i].getY() << ", " <<
1228  nn[i].getZ() << ") on " <<
1229  nn[i].elecCompt().element()->getName() <<
1230  "\n";
1231  }
1232  }
1233  */
1234 
1235  mc->indexToSpace( 72, x, y, z );
1236  near = mc->nearest( x, y, z, index );
1237  assert( index == 72 );
1238  assert( doubleEq( near, 0 ) );
1239 
1240  for( unsigned int i = 0; i < ndc; ++i ) {
1241  mc->indexToSpace( i, x, y, z );
1242  near = mc->nearest( x, y, z + diffLength/10.0, index );
1243  assert( index == i );
1244  assert( near == diffLength / 10.0 );
1245  }
1246 
1248  // Test the area calculation
1249 
1250  CubeMesh cube;
1251  cube.setPreserveNumEntries( 0 );
1252  vector< double > coords( 9, 0.0 );
1253  coords[0] = -10; // X0
1254  coords[1] = -10; // Y0
1255  coords[2] = -10; // Z0
1256 
1257  coords[3] = 10; // X1
1258  coords[4] = 10; // Y1
1259  coords[5] = 10; // Z1
1260  coords[6] = 20; // DX
1261  coords[7] = 20; // DY
1262  coords[8] = 20; // DZ
1263  cube.innerSetCoords( coords );
1264  vector< VoxelJunction > vj;
1265  mc->matchMeshEntries( &cube, vj );
1266  assert( vj.size() == mc->innerGetNumEntries() );
1267  double area = 0.0;
1268  for ( unsigned int i = 0; i < vj.size(); ++i ) {
1269  assert( vj[i].first == i );
1270  assert( vj[i].second == 0 ); // Only one cube entry
1271  area += vj[i].diffScale;
1272  /*
1273  double r = cm.getR0() +
1274  ( cm.getR1() - cm.getR0() ) *
1275  cm.getDiffLength() * ( 0.5 + i ) / cm.getTotLength();
1276  double a = r * cm.getDiffLength() * 2 * PI;
1277  //assert( doubleApprox( vj[i].diffScale, a ) );
1278  // cout << i << ". mesh: " << vj[i].diffScale << ", calc: " << a << endl;
1279  assert( fabs( vj[i].diffScale - a ) < 0.5 );
1280  */
1281  }
1282  double a2 = dia * dia * PI +
1283  4 * len * dia * PI / sqrt( 2.0 ) +
1284  4 * len * dia * PI / 2.0 +
1285  8 * len * dia * PI / (2.0 * sqrt( 2.0 ) );
1286 
1287  // Scaling is needed otherwise doubleApprox treats both as zero.
1288  assert( doubleApprox( area * 1e8, a2 * 1e8 ) );
1289 
1290 
1292 
1293  shell->doDelete( cell );
1294  shell->doDelete( nm );
1295  cout << "." << flush;
1296 }
1297 #endif
1298 
1299 // Assorted definitions from CubeMesh.cpp
1300 static const unsigned int EMPTY = ~0;
1301 static const unsigned int SURFACE = ~1;
1302 static const unsigned int ABUT = ~2;
1303 static const unsigned int MULTI = ~3;
1304 typedef pair< unsigned int, unsigned int > PII;
1305 extern void setIntersectVoxel(
1306  vector< PII >& intersect,
1307  unsigned int ix, unsigned int iy, unsigned int iz,
1308  unsigned int nx, unsigned int ny, unsigned int nz,
1309  unsigned int meshIndex );
1310 
1311 extern void checkAbut(
1312  const vector< PII >& intersect,
1313  unsigned int ix, unsigned int iy, unsigned int iz,
1314  unsigned int nx, unsigned int ny, unsigned int nz,
1315  unsigned int meshIndex,
1316  vector< VoxelJunction >& ret );
1317 
1318 void testIntersectVoxel()
1319 {
1341  unsigned int nx = 5;
1342  unsigned int ny = 3;
1343  unsigned int nz = 1;
1344  vector< PII > intersect( nx * ny * nz, PII(
1345  CubeMesh::EMPTY, CubeMesh::EMPTY ) );
1346  unsigned int meshIndex = 0;
1347  setIntersectVoxel( intersect, 1, 0, 0, nx, ny, nz, meshIndex++ );
1348  setIntersectVoxel( intersect, 2, 0, 0, nx, ny, nz, meshIndex++ );
1349  setIntersectVoxel( intersect, 3, 0, 0, nx, ny, nz, meshIndex++ );
1350  setIntersectVoxel( intersect, 1, 1, 0, nx, ny, nz, meshIndex++ );
1351  setIntersectVoxel( intersect, 1, 2, 0, nx, ny, nz, meshIndex++ );
1352  setIntersectVoxel( intersect, 2, 2, 0, nx, ny, nz, meshIndex++ );
1353  setIntersectVoxel( intersect, 3, 2, 0, nx, ny, nz, meshIndex++ );
1354 
1355  assert( intersect[0].first == 0 &&
1356  intersect[0].second == CubeMesh::ABUTX );
1357  assert( intersect[1].first == 0 &&
1358  intersect[1].second == CubeMesh::SURFACE );
1359  assert( intersect[2].first == 1 &&
1360  intersect[2].second == CubeMesh::SURFACE );
1361  assert( intersect[3].first == 2 &&
1362  intersect[3].second == CubeMesh::SURFACE );
1363  assert( intersect[4].first == 2 &&
1364  intersect[4].second == CubeMesh::ABUTX );
1365 
1366  assert( intersect[5].first == 3 &&
1367  intersect[5].second == CubeMesh::ABUTX );
1368  assert( intersect[6].first == 3 &&
1369  intersect[6].second == CubeMesh::SURFACE );
1370  assert( intersect[7].first == 1 &&
1371  intersect[7].second == CubeMesh::MULTI );
1372  assert( intersect[8].first == 2 &&
1373  intersect[8].second == CubeMesh::MULTI );
1374  assert( intersect[9].first == EMPTY &&
1375  intersect[9].second == CubeMesh::EMPTY );
1376 
1377  assert( intersect[10].first == 4 &&
1378  intersect[10].second == CubeMesh::ABUTX );
1379  assert( intersect[11].first == 4 &&
1380  intersect[11].second == CubeMesh::SURFACE );
1381  assert( intersect[12].first == 5 &&
1382  intersect[12].second == CubeMesh::SURFACE );
1383  assert( intersect[13].first == 6 &&
1384  intersect[13].second == CubeMesh::SURFACE );
1385  assert( intersect[14].first == 6 &&
1386  intersect[14].second == CubeMesh::ABUTX );
1387 
1388  // Next: test out checkAbut.
1389  vector< VoxelJunction > ret;
1390  checkAbut( intersect, 0, 0, 0, nx, ny, nz, 1234, ret );
1391  assert( ret.size() == 1 );
1392  assert( ret[0].first == 0 && ret[0].second == 1234 );
1393  ret.clear();
1394  // The ones below are either SURFACE or EMPTY and should not add points
1395  checkAbut( intersect, 1, 0, 0, nx, ny, nz, 1234, ret );
1396  checkAbut( intersect, 2, 0, 0, nx, ny, nz, 1234, ret );
1397  checkAbut( intersect, 3, 0, 0, nx, ny, nz, 1234, ret );
1398  checkAbut( intersect, 1, 1, 0, nx, ny, nz, 1234, ret );
1399  checkAbut( intersect, 4, 1, 0, nx, ny, nz, 1234, ret );
1400  checkAbut( intersect, 1, 2, 0, nx, ny, nz, 1234, ret );
1401  checkAbut( intersect, 2, 2, 0, nx, ny, nz, 1234, ret );
1402  checkAbut( intersect, 3, 2, 0, nx, ny, nz, 1234, ret );
1403  assert( ret.size() == 0 );
1404  checkAbut( intersect, 2, 1, 0, nx, ny, nz, 9999, ret );
1405  assert( ret.size() == 3 );
1406  assert( ret[0].first == 3 && ret[0].second == 9999 );
1407  assert( ret[1].first == 1 && ret[1].second == 9999 );
1408  assert( ret[2].first == 5 && ret[1].second == 9999 );
1409  ret.clear();
1410  checkAbut( intersect, 3, 1, 0, nx, ny, nz, 8888, ret );
1411  assert( ret.size() == 2 );
1412  assert( ret[0].first == 2 && ret[0].second == 8888 );
1413  assert( ret[1].first == 6 && ret[1].second == 8888 );
1414  ret.clear();
1415  checkAbut( intersect, 4, 0, 0, nx, ny, nz, 7777, ret );
1416  checkAbut( intersect, 0, 1, 0, nx, ny, nz, 6666, ret );
1417  checkAbut( intersect, 0, 2, 0, nx, ny, nz, 5555, ret );
1418  checkAbut( intersect, 4, 2, 0, nx, ny, nz, 4444, ret );
1419  assert( ret.size() == 4 );
1420  assert( ret[0].first == 2 && ret[0].second == 7777 );
1421  assert( ret[1].first == 3 && ret[1].second == 6666 );
1422  assert( ret[2].first == 4 && ret[2].second == 5555 );
1423  assert( ret[3].first == 6 && ret[3].second == 4444 );
1424 
1425  cout << "." << flush;
1426 }
1427 
1428 void testCubeMeshFillTwoDimSurface()
1429 {
1430  CubeMesh cm;
1431  vector< double > coords( 9, 0.0 );
1432  coords[3] = 5.0;
1433  coords[4] = 3.0;
1434  coords[5] = 1.0;
1435  coords[6] = coords[7] = coords[8] = 1.0;
1436  cm.setPreserveNumEntries( false );
1437  cm.innerSetCoords( coords );
1438  assert( cm.numDims() == 2 );
1439  const vector< unsigned int >& surface = cm.surface();
1440  assert( surface.size() == 15 );
1441  for ( unsigned int i = 0; i < 15; ++i ) {
1442  assert( surface[i] == i );
1443  }
1444  cout << "." << flush;
1445 }
1446 
1447 void testCubeMeshFillThreeDimSurface()
1448 {
1449  cout << "." << flush;
1450 }
1451 
1452 void testCubeMeshJunctionTwoDimSurface()
1453 {
1463  CubeMesh cm1;
1464  vector< double > coords( 9, 0.0 );
1465  coords[3] = 5.0;
1466  coords[4] = 3.0;
1467  coords[5] = 1.0;
1468  coords[6] = coords[7] = coords[8] = 1.0;
1469  cm1.setPreserveNumEntries( false );
1470  cm1.innerSetCoords( coords );
1471  vector< unsigned int > surface = cm1.surface();
1472  assert( surface.size() == 15 );
1473 
1474  CubeMesh cm2;
1475  coords[0] = 5.0;
1476  coords[1] = -1.0;
1477  coords[2] = 0.0;
1478  coords[3] = 7.0;
1479  coords[4] = 4.0;
1480  coords[5] = 1.0;
1481  coords[6] = coords[7] = coords[8] = 1.0;
1482  cm2.setPreserveNumEntries( false );
1483  cm2.innerSetCoords( coords );
1484  const vector< unsigned int >& surface2 = cm2.surface();
1485  assert( surface2.size() == 10 );
1486 
1487  vector< VoxelJunction > ret;
1488  cm1.matchCubeMeshEntries( &cm2, ret );
1489  assert( ret.size() == 3 );
1490 
1491  assert( ret[0].first == 4 );
1492  assert( ret[0].second == 2 );
1493  assert( ret[1].first == 9 );
1494  assert( ret[1].second == 4 );
1495  assert( ret[2].first == 14 );
1496  assert( ret[2].second == 6 );
1497 
1510  // Trimming cm1. At this point we don't assume automatic updates of
1511  // the m2s, s2m and surface vectors when any of them is changed.
1512  vector< unsigned int > m2s = cm1.getMeshToSpace();
1513  assert( m2s.size() == 15 );
1514  m2s.resize( 14 );
1515  cm1.setMeshToSpace( m2s );
1516  vector< unsigned int > s2m = cm1.getSpaceToMesh();
1517  assert( s2m.size() == 15 );
1518  s2m[14] = ~0;
1519  cm1.setSpaceToMesh( s2m );
1520  surface.resize( 4 );
1521  // As a shortcut, just assign the places near the junction
1522  // Note that the indices are spaceIndices.
1523  surface[0] = 3;
1524  surface[1] = 4;
1525  surface[2] = 9;
1526  surface[3] = 13;
1527  cm1.setSurface( surface );
1528 
1529  // Trimming cm2.
1530  m2s = cm2.getMeshToSpace();
1531  assert( m2s.size() == 10 );
1532  m2s.resize( 8 );
1533  m2s[0] = 1;
1534  for ( unsigned int i = 1; i < 8; ++i )
1535  m2s[i] = i + 2;
1536  cm2.setMeshToSpace( m2s );
1537  s2m.clear();
1538  s2m.resize( 10, ~0 );
1539  for ( unsigned int i = 0; i < 8; ++i )
1540  s2m[ m2s[i] ] = i;
1541  cm2.setSpaceToMesh( s2m );
1542  // As a shortcut, just assign the places near the junction
1543  // Note that the indices are spaceIndices.
1544  surface[0] = 3;
1545  surface[1] = 4;
1546  surface[2] = 6;
1547  surface[3] = 8;
1548  cm2.setSurface( surface );
1549 
1550  // Now test it out.
1551  ret.resize( 0 );
1552  cm1.matchCubeMeshEntries( &cm2, ret );
1553  assert( ret.size() == 1 );
1554  assert( ret[0].first == 9 );
1555  assert( ret[0].second == 2 );
1556 
1557  cout << "." << flush;
1558 }
1559 
1560 void testCubeMeshJunctionDiffSizeMesh()
1561 {
1574  CubeMesh cm1;
1575  vector< double > coords( 9, 0.0 );
1576  coords[3] = 5.0;
1577  coords[4] = 3.0;
1578  coords[5] = 1.0;
1579  coords[6] = coords[7] = coords[8] = 1.0;
1580  cm1.setPreserveNumEntries( false );
1581  cm1.innerSetCoords( coords );
1582  vector< unsigned int > surface = cm1.surface();
1583  assert( surface.size() == 15 );
1584 
1585  CubeMesh cm2;
1586  coords[0] = 5.0;
1587  coords[1] = -0.5;
1588  coords[2] = 0.0;
1589  coords[3] = 7.0;
1590  coords[4] = 3.5;
1591  coords[5] = 0.5;
1592  coords[6] = 1.0;
1593  coords[7] = 0.5;
1594  coords[8] = 0.5;
1595  cm2.setPreserveNumEntries( false );
1596  cm2.innerSetCoords( coords );
1597  const vector< unsigned int >& surface2 = cm2.surface();
1598  assert( surface2.size() == 16 );
1599 
1600  vector< VoxelJunction > ret;
1601  cm1.matchCubeMeshEntries( &cm2, ret );
1602  assert( ret.size() == 6 );
1603 
1604  assert( ret[0].first == 4 );
1605  assert( ret[0].second == 2 );
1606  assert( ret[1].first == 4 );
1607  assert( ret[1].second == 4 );
1608  assert( ret[2].first == 9 );
1609  assert( ret[2].second == 6 );
1610  assert( ret[3].first == 9 );
1611  assert( ret[3].second == 8 );
1612  assert( ret[4].first == 14 );
1613  assert( ret[4].second == 10 );
1614  assert( ret[5].first == 14 );
1615  assert( ret[5].second == 12 );
1616 
1617  cout << "." << flush;
1618 }
1619 
1620 void testCubeMeshJunctionThreeDimSurface()
1621 {
1622  cout << "." << flush;
1623 }
1624 
1636 void testCubeMeshMultiJunctionTwoD()
1637 {
1638  CubeMesh A;
1639  vector< double > coords( 9, 0.0 );
1640  coords[3] = 10e-6;
1641  coords[4] = 10e-6;
1642  coords[5] = 10e-6;
1643  coords[6] = coords[7] = coords[8] = 10e-6;
1644  A.setPreserveNumEntries( false );
1645  A.innerSetCoords( coords );
1646  vector< unsigned int > surface = A.surface();
1647  assert( surface.size() == 1 );
1648  assert( surface[0] == 0 );
1649 
1650  CubeMesh B;
1651  coords[0] = -30e-6; coords[1] = 0; coords[2] = 0;
1652  coords[3] = 0; coords[4] = 10e-6; coords[5] = 10e-6;
1653  coords[6] = coords[7] = coords[8] = 10e-6;
1654  B.setPreserveNumEntries( false );
1655  B.innerSetCoords( coords );
1656  surface = B.surface();
1657  assert( surface.size() == 3 );
1658  assert( surface[0] == 0 );
1659  assert( surface[1] == 1 );
1660  assert( surface[2] == 2 );
1661 
1662  CubeMesh D;
1663  coords[0] = 0; coords[1] = 10e-6; coords[2] = 0;
1664  coords[3] = 10e-6; coords[4] = 30e-6; coords[5] = 10e-6;
1665  coords[6] = coords[7] = coords[8] = 10e-6;
1666  D.setPreserveNumEntries( false );
1667  D.innerSetCoords( coords );
1668  surface = D.surface();
1669  assert( surface.size() == 2 );
1670  assert( surface[0] == 0 );
1671  assert( surface[1] == 1 );
1672 
1673  CubeMesh C;
1674  coords[0] = -30e-6; coords[1] = -10e-6; coords[2] = 0;
1675  coords[3] = 20e-6; coords[4] = 0; coords[5] = 10e-6;
1676  coords[6] = coords[7] = coords[8] = 10e-6;
1677  C.setPreserveNumEntries( false );
1678  C.innerSetCoords( coords );
1679  surface = C.surface();
1680  assert( surface.size() == 5 );
1681  assert( surface[0] == 0 );
1682  assert( surface[1] == 1 );
1683  assert( surface[2] == 2 );
1684  assert( surface[3] == 3 );
1685  assert( surface[4] == 4 );
1686 
1687 
1688  cout << "." << flush;
1689 }
1690 
1691 void testVec()
1692 {
1693  Vec i( 1, 0, 0 );
1694  Vec j( 0, 1, 0 );
1695  Vec k( 0, 0, 1 );
1696 
1697  assert( doubleEq( i.dotProduct( j ), 0.0 ) );
1698  assert( doubleEq( j.dotProduct( k ), 0.0 ) );
1699  assert( doubleEq( j.dotProduct( j ), 1.0 ) );
1700 
1701  assert( i.crossProduct( j ) == k );
1702 
1703  Vec u, v;
1704  i.orthogonalAxes( u, v );
1705  assert( u == j );
1706  assert( v == k );
1707 
1708  cout << "." << flush;
1709 }
1710 
1711 #if 0
1712 void testSpineEntry()
1713 {
1714  Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
1715  // Build a cell
1716  Id cell = shell->doCreate( "Neutral", Id(), "cell", 1 );
1717 // Id makeCompt( Id parentCompt, Id parentObj, string name, double len, double dia, double theta )
1718  Id neck = makeCompt( Id(), cell, "neck", 1e-6, 1e-7, 0 );
1719  Id head = makeCompt( neck, cell, "head", 1e-6, 1e-6, 0 );
1720 
1721  SpineEntry se( neck, head, 1234 );
1722 
1723  assert( se.parent() == 1234 );
1724  assert( se.shaftId() == neck );
1725  assert( se.headId() == head );
1726  assert( doubleEq( se.volume(), 1e-18 * PI/4.0 ) );
1727  double x, y, z;
1728  se.mid( x, y, z );
1729  assert( doubleEq( x, 1.5e-6 ) );
1730  assert( doubleEq( y, 0.0 ) );
1731  assert( doubleEq( z, 0.0 ) );
1732 
1734  // Test matchCubeMeshEntries.
1735  vector< VoxelJunction > vj;
1736  CubeMesh cube;
1737  cube.setPreserveNumEntries( 0 );
1738  vector< double > coords( 9, 0.0 );
1739  coords[0] = -10; // X0
1740  coords[1] = -10; // Y0
1741  coords[2] = -10; // Z0
1742 
1743  coords[3] = 10; // X1
1744  coords[4] = 10; // Y1
1745  coords[5] = 10; // Z1
1746  coords[6] = 20; // DX
1747  coords[7] = 20; // DY
1748  coords[8] = 20; // DZ
1749  cube.innerSetCoords( coords );
1750  se.matchCubeMeshEntriesToHead( &cube, 1234, 0.1, vj );
1751  assert( vj.size() == 1 );
1752  assert( vj[0].first == 1234 );
1753  assert( vj[0].second == 0 );
1754  assert( doubleApprox( vj[0].diffScale * 1e10, 1e-12 * PI * 1e10 ) );
1755 
1756  vj.clear();
1757  se.matchCubeMeshEntriesToPSD( &cube, 4321, 0.1, vj );
1758  assert( vj.size() == 1 );
1759  assert( vj[0].first == 4321 );
1760  assert( vj[0].second == 0 );
1761  assert( doubleApprox( vj[0].diffScale * 1e10, 1e-12 * PI * 0.25 * 1e10 ) );
1762 
1764  shell->doDelete( cell );
1765  cout << "." << flush;
1766 }
1767 
1768 void testSpineAndPsdMesh()
1769 {
1770  Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
1771  // Build a linear cylindrical cell, no tapering.
1772  Id cell = shell->doCreate( "Neutral", Id(), "cell", 1 );
1773  unsigned int numCompts = 500;
1774  double dia = 20e-6; // metres
1775  double diffLength = 0.5e-6; // metres
1776  double len = diffLength * numCompts;
1777  unsigned int numSpines = 10;
1778  Id soma = makeCompt( Id(), cell, "soma", dia, dia, 90 );
1779  Field< double >::set( soma, "y0", -dia );
1780  Field< double >::set( soma, "y", 0 );
1781  Id dend = makeCompt( soma, cell, "dend", len, dia, 0 );
1782 
1783  for ( unsigned int i = 0; i < numSpines; ++i ) {
1784  double frac = i / static_cast< double >( numSpines );
1785  makeSpine( dend, cell, i, frac, 1.0e-6, 1.0e-6, i * 30.0 );
1786  }
1787 
1788  Id nm = shell->doCreate( "NeuroMesh", Id(), "neuromesh", 1 );
1789  Field< bool >::set( nm, "separateSpines", true );
1790  Field< double >::set( nm, "diffLength", diffLength );
1791  Field< string >::set( nm, "geometryPolicy", "cylinder" );
1792  Id sm = shell->doCreate( "SpineMesh", Id(), "spinemesh", 1 );
1793  ObjId mid = shell->doAddMsg(
1794  "OneToOne", nm, "spineListOut", sm, "spineList" );
1795  assert( !mid.bad() );
1796  Id pm = shell->doCreate( "PsdMesh", Id(), "psdmesh", 1 );
1797  mid = shell->doAddMsg(
1798  "OneToOne", nm, "psdListOut", pm, "psdList" );
1799  assert( !mid.bad() );
1800  Field< Id >::set( nm, "cell", cell );
1801  unsigned int ns = Field< unsigned int >::get( nm, "numSegments" );
1802  assert( ns == 2 ); // soma and dend
1803  unsigned int ndc = Field< unsigned int >::get( nm, "numDiffCompts" );
1804  assert( ndc == numCompts + floor( dia / diffLength + 0.5 ) );
1805 
1806  unsigned int sdc = Field< unsigned int >::get( sm, "num_mesh" );
1807  assert( sdc == numSpines );
1808 
1809  // Should now iterate through the spines to find their 'nearest' on
1810  // the dendrite.
1811  SpineMesh* s = reinterpret_cast< SpineMesh* >( sm.eref().data() );
1812  for ( unsigned int i = 0; i < sdc; ++i ) {
1813  assert( s->spines()[i].parent() == 40 + i * numCompts/numSpines );
1814  assert( doubleEq( s->spines()[i].volume(), PI * 0.25e-12 ) );
1815  double x = i * diffLength * numCompts/numSpines;
1816  unsigned int index = 0;
1817  s->nearest( x, 0, 0, index );
1818  assert( index == i );
1819  }
1820 
1821  // Check coupling to Neuromesh.
1822  vector< VoxelJunction > ret;
1823  NeuroMesh* nmesh = reinterpret_cast< NeuroMesh* >( nm.eref().data() );
1824  s->matchMeshEntries( nmesh, ret );
1825  assert( ret.size() == numSpines );
1826  for ( unsigned int i = 0; i < numSpines; ++i ) {
1827  assert( ret[i].first == i );
1828  assert( ret[i].second == s->spines()[i].parent() );
1829  assert( doubleEq( ret[i].diffScale, PI * 0.25e-14 / 1.5e-6 ) );
1830  }
1831 
1832  // Check coupling to CubeMesh.
1833  CubeMesh cube;
1834  cube.setPreserveNumEntries( 0 );
1835  vector< double > coords( 9, 0.0 );
1836  coords[0] = -len/2.0 - 1.2e-5; // X0
1837  coords[1] = -1e-5; // Y0
1838  coords[2] = -1e-5; // Z0
1839 
1840  coords[3] = len * 2.0 - 1.2e-5; // X1
1841  coords[4] = 1.5e-5; // Y1
1842  coords[5] = 1.5e-5; // Z1
1843  coords[6] = 2.5e-5; // DX
1844  coords[7] = 2.5e-5; // DY
1845  coords[8] = 2.5e-5; // DZ
1846  cube.innerSetCoords( coords );
1847  ret.clear();
1848  s->matchMeshEntries( &cube, ret );
1849  assert( ret.size() == numSpines );
1850  for ( unsigned int i = 0; i < ret.size(); ++i ) {
1851  assert( ret[i].first == i );
1852  unsigned int cubeIndex = i + 5;
1853  assert( ret[i].second == cubeIndex );
1854  assert( doubleApprox( ret[i].diffScale * 1e10, PI * 1e-12 * 1e10 ) );
1855  /*
1856  cout << i << " cubeIndex=" << cubeIndex << ", (" <<
1857  ret[i].first << ", " << ret[i].second << ") : " <<
1858  ret[i].diffScale << "\n";
1859  */
1860  }
1862  // Check PSD mesh
1863  unsigned int pdc = Field< unsigned int >::get( pm, "num_mesh" );
1864  Id pe( pm.value() + 1 );
1865  assert( pdc == numSpines );
1866  PsdMesh* p = reinterpret_cast< PsdMesh* >( pm.eref().data() );
1867  for ( unsigned int i = 0; i < pdc; ++i ) {
1868  assert( p->parent( i ) == i );
1869  ObjId oi( pe, i );
1870  double area = Field< double >::get( oi, "volume" );
1871  assert( doubleEq( area, 1e-12 * 0.25 * PI ) );
1872  unsigned int dim = Field< unsigned int >::get( oi, "dimensions" );
1873  assert( dim == 2 );
1874  double x = i * diffLength * numCompts/numSpines;
1875  unsigned int index = 0;
1876  s->nearest( x, 0, 0, index );
1877  assert( index == i );
1878  }
1879  ret.clear();
1880  p->matchMeshEntries( s, ret );
1881  assert( ret.size() == numSpines );
1882  for ( unsigned int i = 0; i < ret.size(); ++i ) {
1883  assert( ret[i].first == i );
1884  assert( ret[i].second == i );
1885  assert( doubleApprox( ret[i].diffScale * 1e10,
1886  1e10 * 0.25 * PI * 1e-12 / 0.5e-6 ) );
1887  }
1888  ret.clear();
1889  p->matchMeshEntries( &cube, ret );
1890  assert( ret.size() == numSpines );
1891  for ( unsigned int i = 0; i < ret.size(); ++i ) {
1892  assert( ret[i].first == i );
1893  unsigned int cubeIndex = i + 5;
1894  assert( ret[i].second == cubeIndex );
1895  assert( doubleApprox( ret[i].diffScale * 1e10, 0.25 * PI * 1e-12 * 1e10 ) );
1896  }
1897 
1898 
1899  shell->doDelete( cell );
1900  shell->doDelete( nm );
1901  shell->doDelete( sm );
1902  shell->doDelete( pm );
1903  cout << "." << flush;
1904 }
1905 
1906 #include "../shell/Wildcard.h"
1907 void testNeuroNodeTree()
1908 {
1909  Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
1910  // Build a cell
1911  Id cell = shell->doCreate( "Neutral", Id(), "cell", 1 );
1912  double len = 10e-6; // metres
1913  double dia = 1e-6; // metres
1914 
1915  buildBranchingCell( cell, len, dia );
1916 
1917  vector< NeuroNode > nodes;
1918  vector< Id > elist;
1919  wildcardFind( "/cell/##", elist );
1920  NeuroNode::buildTree( nodes, elist );
1921  assert( nodes.size() == 17 );
1922  unsigned int i = 0;
1923  // Here we want the function to generate a depth-first tree.
1924  assert( nodes[i++].elecCompt().element()->getName() == "soma" );
1925  assert( nodes[i++].elecCompt().element()->getName() == "d1" );
1926  assert( nodes[i++].elecCompt().element()->getName() == "d1a" );
1927  assert( nodes[i++].elecCompt().element()->getName() == "d11" );
1928  assert( nodes[i++].elecCompt().element()->getName() == "d111" );
1929  assert( nodes[i++].elecCompt().element()->getName() == "d112" );
1930  assert( nodes[i++].elecCompt().element()->getName() == "d12" );
1931  assert( nodes[i++].elecCompt().element()->getName() == "d121" );
1932  assert( nodes[i++].elecCompt().element()->getName() == "d122" );
1933 
1934  assert( nodes[i++].elecCompt().element()->getName() == "d2" );
1935  assert( nodes[i++].elecCompt().element()->getName() == "d2a" );
1936  assert( nodes[i++].elecCompt().element()->getName() == "d21" );
1937  assert( nodes[i++].elecCompt().element()->getName() == "d211" );
1938  assert( nodes[i++].elecCompt().element()->getName() == "d212" );
1939  assert( nodes[i++].elecCompt().element()->getName() == "d22" );
1940  assert( nodes[i++].elecCompt().element()->getName() == "d221" );
1941  assert( nodes[i++].elecCompt().element()->getName() == "d222" );
1942 
1943  /*
1944  Id soma = makeCompt( Id(), cell, "soma", dia, dia, 90 );
1945  Id dend = makeCompt( Id(), cell, "dend", len, dia, 0 );
1946  for ( unsigned int i = 0; i < numSpines; ++i ) {
1947  double frac = i / static_cast< double >( numSpines );
1948  makeSpine( dend, cell, i, frac, 1.0e-6, 1.0e-6, i * 30.0 );
1949  }
1950  */
1952  // Now to make a totally random cell, with a complete mash of
1953  // messaging, having circular linkage even.
1955  shell->doDelete( cell );
1956  cell = shell->doCreate( "Neutral", Id(), "cell", 1 );
1957  vector< Id > compt( 10 );
1958  for ( unsigned int i = 0; i < compt.size(); ++i ) {
1959  stringstream ss;
1960  ss << "compt" << i;
1961  string name = ss.str();
1962  compt[i] = shell->doCreate( "SymCompartment", cell, name, 1 );
1963  double dia = 50.0 - (i - 3) * (i - 3 );
1964  Field< double >::set( compt[i], "diameter", dia );
1965  }
1966  shell->doAddMsg( "Single", compt[0], "proximal", compt[1], "distal" );
1967  shell->doAddMsg( "Single", compt[1], "proximal", compt[2], "distal" );
1968  shell->doAddMsg( "Single", compt[2], "distal", compt[3], "proximal" );
1969  shell->doAddMsg( "Single", compt[3], "sibling", compt[2], "sibling" );
1970  shell->doAddMsg( "Single", compt[3], "sphere", compt[4], "proximalOnly" );
1971  shell->doAddMsg( "Single", compt[4], "distal", compt[5], "proximal" );
1972  shell->doAddMsg( "Single", compt[4], "distal", compt[6], "proximal" );
1973  shell->doAddMsg( "Single", compt[4], "axial", compt[7], "raxial" );
1974  shell->doAddMsg( "Single", compt[7], "cylinder", compt[8], "proximalOnly" );
1975  shell->doAddMsg( "Single", compt[7], "cylinder", compt[9], "proximalOnly" );
1976  shell->doAddMsg( "Single", compt[9], "proximal", compt[0], "distal" );
1977 
1978  nodes.clear();
1979  elist.clear();
1980  wildcardFind( "/cell/##", elist );
1981  NeuroNode::buildTree( nodes, elist );
1982  assert( nodes.size() == 10 );
1983  i = 0;
1984  // The fallback for the soma is the compt with biggest diameter.
1985  assert( nodes[i++].elecCompt().element()->getName() == "compt3" );
1986  assert( nodes[i++].elecCompt().element()->getName() == "compt2" );
1987  assert( nodes[i++].elecCompt().element()->getName() == "compt1" );
1988  assert( nodes[i++].elecCompt().element()->getName() == "compt0" );
1989  assert( nodes[i++].elecCompt().element()->getName() == "compt9" );
1990  assert( nodes[i++].elecCompt().element()->getName() == "compt7" );
1991  assert( nodes[i++].elecCompt().element()->getName() == "compt8" );
1992  assert( nodes[i++].elecCompt().element()->getName() == "compt4" );
1993  assert( nodes[i++].elecCompt().element()->getName() == "compt5" );
1994  assert( nodes[i++].elecCompt().element()->getName() == "compt6" );
1995 
1996  // Now I change the biggest compartment to compartment 5.
1997  Field< double >::set( compt[5], "diameter", 51.0 );
1998  nodes.clear();
1999  elist.clear();
2000  wildcardFind( "/cell/##", elist );
2001  NeuroNode::buildTree( nodes, elist );
2002  assert( nodes.size() == 10 );
2003  i = 0;
2004  assert( nodes[i++].elecCompt().element()->getName() == "compt5" );
2005  assert( nodes[i++].elecCompt().element()->getName() == "compt4" );
2006  assert( nodes[i++].elecCompt().element()->getName() == "compt3" );
2007  assert( nodes[i++].elecCompt().element()->getName() == "compt2" );
2008  assert( nodes[i++].elecCompt().element()->getName() == "compt1" );
2009  assert( nodes[i++].elecCompt().element()->getName() == "compt0" );
2010  assert( nodes[i++].elecCompt().element()->getName() == "compt9" );
2011  assert( nodes[i++].elecCompt().element()->getName() == "compt7" );
2012  assert( nodes[i++].elecCompt().element()->getName() == "compt8" );
2013  assert( nodes[i++].elecCompt().element()->getName() == "compt6" );
2015  // Now to make a dendrite with lots of spines and filter them off.
2017  shell->doDelete( cell );
2018  cell = shell->doCreate( "Neutral", Id(), "cell", 1 );
2019  vector< Id > dend( 10 );
2020  vector< Id > neck( 10 );
2021  vector< Id > head( 10 );
2022  for ( unsigned int i = 0; i < neck.size(); ++i ) {
2023  stringstream ss;
2024  ss << "dend" << i;
2025  string name = ss.str();
2026  dend[i] = shell->doCreate( "SymCompartment", cell, name, 1 );
2027  Field< double >::set( dend[i], "diameter", 2.0e-6 );
2028  stringstream ss1;
2029  ss1 << "neck" << i;
2030  name = ss1.str();
2031  neck[i] = shell->doCreate( "SymCompartment", cell, name, 1 );
2032  Field< double >::set( neck[i], "diameter", 0.2e-6 );
2033  stringstream ss2;
2034  ss2 << "head" << i;
2035  name = ss2.str();
2036  head[i] = shell->doCreate( "SymCompartment", cell, name, 1 );
2037  Field< double >::set( head[i], "diameter", 1.0e-6 );
2038  }
2039  for ( unsigned int i = 0; i < neck.size(); ++i ) {
2040  // Just to make things interesting, add the messages staggered.
2041  if ( i > 0 )
2042  shell->doAddMsg( "Single", dend[i], "proximal",
2043  dend[ (i-1)], "distal" );
2044  shell->doAddMsg( "Single", neck[i], "proximalOnly",
2045  dend[ (i+7) % head.size()], "cylinder" );
2046  shell->doAddMsg( "Single", neck[i], "distal",
2047  head[ (i+4) % neck.size() ], "proximal" );
2048  }
2049  nodes.clear();
2050  elist.clear();
2051  wildcardFind( "/cell/##", elist );
2052  NeuroNode::buildTree( nodes, elist );
2053  assert( nodes.size() == 30 );
2054  vector< NeuroNode > orig = nodes;
2055  vector< Id > shaftId;
2056  vector< Id > headId;
2057  vector< unsigned int > parent;
2058  NeuroNode::filterSpines( nodes, shaftId, headId, parent );
2059  assert( nodes.size() == 10 );
2060  /*
2061  cout << endl;
2062  for ( unsigned int i = 0; i < 10; ++i ) {
2063  cout << "nodes[" << i << "].pa = " << nodes[i].parent() <<
2064  ", Id = " << nodes[i].elecCompt() << endl;
2065  }
2066  for ( unsigned int i = 0; i < 30; ++i ) {
2067  cout << "orig[" << i << "].pa = " << orig[i].parent() <<
2068  ", Id = " << orig[i].elecCompt() << endl;
2069  }
2070  */
2071  // Check that the correct dend is parent of each shaft, and
2072  // correct shaft is parent of each head.
2073  for ( unsigned int i = 0; i < parent.size(); ++i ) {
2074  unsigned int j = (shaftId[i].value() - neck[0].value())/3;
2075  assert( parent[i] == (j+7) % 10 );
2076  unsigned int k = (headId[i].value() - head[0].value())/3;
2077  assert( k == (j+4) % 10 );
2078  }
2079 
2080  shell->doDelete( cell );
2081  cout << "." << flush;
2082 }
2083 
2084 void testCellPortion()
2085 {
2086  Shell* shell = reinterpret_cast< Shell* >( Id().eref().data() );
2087  // Build a linear cylindrical cell, no tapering.
2088  Id cell = shell->doCreate( "Neutral", Id(), "cell", 1 );
2089  unsigned int numCompts = 50;
2090  double dia = 20e-6; // metres
2091  double diffLength = 0.5e-6; // metres
2092  double len = diffLength * numCompts;
2093  unsigned int numSpines = 10;
2094  unsigned int numSpines2 = 25;
2095  Id soma = makeCompt( Id(), cell, "soma", dia, dia, 90 );
2096  Field< double >::set( soma, "y0", -dia );
2097  Field< double >::set( soma, "y", 0 );
2098  Id dend0 = makeCompt( soma, cell, "dend0", len, dia, 0 );
2099  Id dend1 = makeCompt( dend0, cell, "dend1", len, dia, 0 );
2100  Id dend2 = makeCompt( dend1, cell, "dend2", len, dia, 0 );
2101 
2102  for ( unsigned int i = 0; i < numSpines; ++i ) {
2103  double frac = i / static_cast< double >( numSpines );
2104  makeSpine( dend1, cell, i, frac, 1.0e-6, 1.0e-6, i * 30.0 );
2105  }
2106  for ( unsigned int i = 0; i < numSpines2; ++i ) {
2107  double frac = i / static_cast< double >( numSpines2 );
2108  makeSpine( dend2, cell, i + numSpines, frac, 1.0e-6, 1.0e-6, i * 30.0 );
2109  }
2110 
2111  Id nm = shell->doCreate( "NeuroMesh", Id(), "neuromesh", 1 );
2112  Field< bool >::set( nm, "separateSpines", true );
2113  Field< double >::set( nm, "diffLength", diffLength );
2114  Field< string >::set( nm, "geometryPolicy", "cylinder" );
2115  Id sm = shell->doCreate( "SpineMesh", Id(), "spinemesh", 1 );
2116  ObjId mid = shell->doAddMsg(
2117  "OneToOne", nm, "spineListOut", sm, "spineList" );
2118  assert( !mid.bad() );
2119  Id pm = shell->doCreate( "PsdMesh", Id(), "psdmesh", 1 );
2120  mid = shell->doAddMsg(
2121  "OneToOne", nm, "psdListOut", pm, "psdList" );
2122  assert( !mid.bad() );
2123  // SetGet2< Id, string >::set( nm, "cellPortion", cell, "/cell/dend2,/cell/shaft#,/cell/head#" );
2124  SetGet2< Id, string >::set( nm, "cellPortion", cell, "/cell/dend1,/cell/dend2,/cell/shaft1?,/cell/head1?,/cell/shaft3,cell/head3" );
2125  unsigned int ns = Field< unsigned int >::get( nm, "numSegments" );
2126  assert( ns == 2 ); // dend1 + dend2 only
2127  unsigned int ndc = Field< unsigned int >::get( nm, "numDiffCompts" );
2128  assert( ndc == numCompts * 2 );
2129 
2130  unsigned int sdc = Field< unsigned int >::get( sm, "num_mesh" );
2131  assert( sdc == 11 ); // I've selected only those in the teens plus #3
2132 
2133  shell->doDelete( cell );
2134  cout << "." << flush;
2135 }
2136 #endif
2137 
2138 void testMesh()
2139 {
2140  testVec();
2141  testVolScaling();
2142  // testCylBase();
2143  // testNeuroNode();
2144  // testNeuroNodeTree();
2145  // testCylMesh();
2146  // testMidLevelCylMesh();
2147  testCubeMesh();
2148  testCubeMeshExtendStencil();
2149  // testReMesh(); // Waiting to have pool subdivision propagate.
2150  // testNeuroMeshLinear();
2151  // testNeuroMeshBranching();
2152  // testIntersectVoxel();
2153  testCubeMeshFillTwoDimSurface();
2154  testCubeMeshFillThreeDimSurface();
2155  testCubeMeshJunctionTwoDimSurface();
2156  testCubeMeshJunctionThreeDimSurface();
2157  testCubeMeshJunctionDiffSizeMesh();
2158  testCubeMeshMultiJunctionTwoD();
2159  // testSpineEntry();
2160  // testSpineAndPsdMesh();
2161  // testCellPortion();
2162 }
testReMesh
void testReMesh()
Definition: testMesh.cpp:792
wildcardFind
int wildcardFind(const string &path, vector< ObjId > &ret)
Definition: Wildcard.cpp:169
NeuroMesh::getMeshEntryVolume
double getMeshEntryVolume(unsigned int fid) const
Virtual function to return volume of mesh Entry.
Definition: NeuroMesh.cpp:906
CylMesh::getDiffusionArea
vector< double > getDiffusionArea(unsigned int fid) const
Virtual function to return diffusion X-section area.
Definition: CylMesh.cpp:520
Eref::data
char * data() const
Definition: Eref.cpp:41
CylMesh::nearest
double nearest(double x, double y, double z, unsigned int &index) const
Definition: CylMesh.cpp:1075
MeshCompt::getStencilRow
unsigned int getStencilRow(unsigned int meshIndex, const double **entry, const unsigned int **colIndex) const
Definition: MeshCompt.cpp:54
CubeMesh::getNx
unsigned int getNx() const
Definition: CubeMesh.cpp:542
ChemCompt::nearest
virtual double nearest(double x, double y, double z, unsigned int &index) const =0
CubeMesh::getNz
unsigned int getNz() const
Definition: CubeMesh.cpp:565
CubeMesh::getDz
double getDz() const
Definition: CubeMesh.cpp:531
makeCompt
Id makeCompt(Id parentCompt, Id parentObj, string name, double len, double dia, double theta)
Definition: testMesh.cpp:855
NA
const double NA
Definition: consts.cpp:15
ABUT
static const unsigned int ABUT
Definition: testMesh.cpp:1302
CylMesh::getY0
double getY0(const Eref &e) const
Definition: CylMesh.cpp:259
PsdMesh::parent
unsigned int parent(unsigned int index) const
Definition: PsdMesh.cpp:395
CylMesh::getCoordinates
vector< double > getCoordinates(unsigned int fid) const
Virtual function to return coords of mesh Entry.
Definition: CylMesh.cpp:489
CylMesh::getMeshType
unsigned int getMeshType(unsigned int fid) const
Virtual function to return MeshType of specified entry.
Definition: CylMesh.cpp:421
MeshCompt::getNeighbors
vector< unsigned int > getNeighbors(unsigned int fid) const
Looks up stencil to return vector of indices of coupled voxels.
Definition: MeshCompt.cpp:69
ObjId::bad
bool bad() const
Definition: ObjId.cpp:18
Id::value
unsigned int value() const
Definition: Id.cpp:197
PII
pair< unsigned int, unsigned int > PII
Definition: testMesh.cpp:1304
setIntersectVoxel
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: CubeMesh.cpp:1304
NeuroNode::filterSpines
static void filterSpines(vector< NeuroNode > &nodes, vector< Id > &shaftId, vector< Id > &headId, vector< unsigned int > &parent)
Definition: NeuroNode.cpp:620
Field
Definition: SetGet.h:236
CubeMesh::getX0
double getX0() const
Definition: CubeMesh.cpp:441
NeuroNode::buildTree
static void buildTree(vector< NeuroNode > &nodes, vector< ObjId > elist)
Definition: NeuroNode.cpp:558
testVolScaling
void testVolScaling()
Definition: testMesh.cpp:34
CubeMesh::setX1
void setX1(double v)
Definition: CubeMesh.cpp:468
testCubeMeshFillThreeDimSurface
void testCubeMeshFillThreeDimSurface()
Definition: testMesh.cpp:1447
CubeMesh::getCoords
vector< double > getCoords(const Eref &e) const
Definition: CubeMesh.cpp:636
CylMesh::setDiffLength
void setDiffLength(const Eref &e, double v)
Definition: CylMesh.cpp:392
CubeMesh::getX1
double getX1() const
Definition: CubeMesh.cpp:474
makeReacTest
Id makeReacTest()
Definition: testKsolve.cpp:34
CubeMesh::getDx
double getDx() const
Definition: CubeMesh.cpp:507
EMPTY
static const unsigned int EMPTY
Definition: testMesh.cpp:1300
ObjId
Definition: ObjId.h:20
CylMesh::indexToSpace
void indexToSpace(unsigned int index, double &x, double &y, double &z) const
Definition: CylMesh.cpp:1026
CubeMesh::setY0
void setY0(double v)
Definition: CubeMesh.cpp:446
SpineMesh::matchMeshEntries
void matchMeshEntries(const ChemCompt *other, vector< VoxelJunction > &ret) const
Definition: SpineMesh.cpp:486
CubeMesh::setZ0
void setZ0(double v)
Definition: CubeMesh.cpp:457
numEntries
const int numEntries
Definition: proc.cpp:60
Id::eref
Eref eref() const
Definition: Id.cpp:125
CylMesh::setY1
void setY1(const Eref &e, double v)
Definition: CylMesh.cpp:304
CylMesh::getCoords
vector< double > getCoords(const Eref &e) const
Definition: CylMesh.cpp:372
CubeMesh::matchCubeMeshEntries
void matchCubeMeshEntries(const CubeMesh *other, vector< VoxelJunction > &ret) const
Definition: CubeMesh.cpp:1453
CubeMesh::getY1
double getY1() const
Definition: CubeMesh.cpp:485
SpineMesh::spines
const vector< SpineEntry > & spines() const
Definition: SpineMesh.cpp:467
Field::set
static bool set(const ObjId &dest, const string &field, A arg)
Definition: SetGet.h:245
CubeMesh::setZ1
void setZ1(double v)
Definition: CubeMesh.cpp:490
CubeMesh::getNy
unsigned int getNy() const
Definition: CubeMesh.cpp:554
CubeMesh::innerSetCoords
void innerSetCoords(const vector< double > &v)
Definition: CubeMesh.cpp:602
Shell::doCreate
Id doCreate(string type, ObjId parent, string name, unsigned int numData, NodePolicy nodePolicy=MooseBlockBalance, unsigned int preferredNode=1)
Definition: Shell.cpp:181
CubeMesh::getZ0
double getZ0() const
Definition: CubeMesh.cpp:463
Vec::crossProduct
Vec crossProduct(const Vec &other) const
Definition: Vec.cpp:26
CubeMesh::surface
const vector< unsigned int > & surface() const
Utility and test function to read surface.
Definition: CubeMesh.cpp:1111
CubeMesh::MULTI
static const unsigned int MULTI
Definition: CubeMesh.h:290
testCubeMeshJunctionTwoDimSurface
void testCubeMeshJunctionTwoDimSurface()
Definition: testMesh.cpp:1452
testCubeMeshJunctionThreeDimSurface
void testCubeMeshJunctionThreeDimSurface()
Definition: testMesh.cpp:1620
CubeMesh::setSpaceToMesh
void setSpaceToMesh(vector< unsigned int > v)
Definition: CubeMesh.cpp:666
ChemCompt::getNumEntries
unsigned int getNumEntries() const
Definition: ChemCompt.cpp:390
doubleApprox
bool doubleApprox(double x, double y)
Definition: doubleEq.cpp:24
CylMesh::getTotLength
double getTotLength() const
Definition: CylMesh.cpp:406
ChemCompt::getDimensions
unsigned int getDimensions() const
Definition: ChemCompt.cpp:340
CylMesh::innerGetNumEntries
unsigned int innerGetNumEntries() const
Definition: CylMesh.cpp:615
CylMesh::getX1
double getX1(const Eref &e) const
Definition: CylMesh.cpp:299
CylMesh::getMeshDimensions
unsigned int getMeshDimensions(unsigned int fid) const
Virtual function to return dimensions of specified entry.
Definition: CylMesh.cpp:430
doubleEq
bool doubleEq(double x, double y)
Definition: doubleEq.cpp:16
Vec
Definition: Vec.h:13
CylMesh::setZ0
void setZ0(const Eref &e, double v)
Definition: CylMesh.cpp:264
CylMesh::getMeshEntryVolume
double getMeshEntryVolume(unsigned int fid) const
Virtual function to return volume of mesh Entry.
Definition: CylMesh.cpp:477
CubeMesh::getMeshToSpace
vector< unsigned int > getMeshToSpace() const
Definition: CubeMesh.cpp:661
CylMesh::setX0
void setX0(const Eref &e, double v)
Definition: CylMesh.cpp:238
CubeMesh::getY0
double getY0() const
Definition: CubeMesh.cpp:452
CYL
Definition: MeshEntry.h:20
CylMesh::matchMeshEntries
void matchMeshEntries(const ChemCompt *other, vector< VoxelJunction > &ret) const
Definition: CylMesh.cpp:827
testVec
void testVec()
Definition: testMesh.cpp:1691
CylMesh::setX1
void setX1(const Eref &e, double v)
Definition: CylMesh.cpp:291
CubeMesh::setMeshToSpace
void setMeshToSpace(vector< unsigned int > v)
Definition: CubeMesh.cpp:655
CylMesh::getR1
double getR1(const Eref &e) const
Definition: CylMesh.cpp:338
CubeMesh::getSpaceToMesh
vector< unsigned int > getSpaceToMesh() const
Definition: CubeMesh.cpp:672
ChemCompt::innerGetNumEntries
virtual unsigned int innerGetNumEntries() const =0
CylMesh::setR0
void setR0(const Eref &e, double v)
Definition: CylMesh.cpp:277
CubeMesh::getMeshType
unsigned int getMeshType(unsigned int fid) const
Virtual function to return MeshType of specified entry.
Definition: CubeMesh.cpp:772
ChemCompt::indexToSpace
virtual void indexToSpace(unsigned int index, double &x, double &y, double &z) const =0
CubeMesh::SURFACE
static const unsigned int SURFACE
Definition: CubeMesh.h:286
MULTI
static const unsigned int MULTI
Definition: testMesh.cpp:1303
checkAbut
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)
Definition: CubeMesh.cpp:1348
CylMesh::getR0
double getR0(const Eref &e) const
Definition: CylMesh.cpp:285
SURFACE
static const unsigned int SURFACE
Definition: testMesh.cpp:1301
CubeMesh::setY1
void setY1(double v)
Definition: CubeMesh.cpp:479
CylMesh::getZ0
double getZ0(const Eref &e) const
Definition: CylMesh.cpp:272
CubeMesh::getZ1
double getZ1() const
Definition: CubeMesh.cpp:496
Vec::dotProduct
double dotProduct(const Vec &other) const
Definition: Vec.cpp:22
testMesh
void testMesh()
Definition: testMesh.cpp:2138
PsdMesh::matchMeshEntries
void matchMeshEntries(const ChemCompt *other, vector< VoxelJunction > &ret) const
Definition: PsdMesh.cpp:443
MeshCompt::extendStencil
void extendStencil(const ChemCompt *other, const vector< VoxelJunction > &vj)
Add boundary voxels to stencil for cross-solver junctions.
Definition: MeshCompt.cpp:121
SpineMesh::nearest
double nearest(double x, double y, double z, unsigned int &index) const
Definition: SpineMesh.cpp:516
CylMesh::getY1
double getY1(const Eref &e) const
Definition: CylMesh.cpp:312
testIntersectVoxel
void testIntersectVoxel()
Definition: testMesh.cpp:1318
name
static char name[]
Definition: mfield.cpp:401
CylMesh::setR1
void setR1(const Eref &e, double v)
Definition: CylMesh.cpp:330
CylMesh::matchCubeMeshEntries
void matchCubeMeshEntries(const CubeMesh *other, vector< VoxelJunction > &ret) const
Definition: CylMesh.cpp:983
testCubeMeshFillTwoDimSurface
void testCubeMeshFillTwoDimSurface()
Definition: testMesh.cpp:1428
CylMesh::getX0
double getX0(const Eref &e) const
Definition: CylMesh.cpp:246
Shell::doAddMsg
ObjId doAddMsg(const string &msgType, ObjId src, const string &srcField, ObjId dest, const string &destField)
Definition: Shell.cpp:269
testCubeMeshJunctionDiffSizeMesh
void testCubeMeshJunctionDiffSizeMesh()
Definition: testMesh.cpp:1560
CubeMesh::EMPTY
static const unsigned int EMPTY
Definition: CubeMesh.h:285
PI
const double PI
Definition: consts.cpp:12
Shell::doDelete
bool doDelete(ObjId oid)
Definition: Shell.cpp:259
Id
Definition: Id.h:17
ChemCompt::matchMeshEntries
virtual void matchMeshEntries(const ChemCompt *other, vector< VoxelJunction > &ret) const =0
Vec::orthogonalAxes
void orthogonalAxes(Vec &u, Vec &v) const
Generates vectors u and v to form a mutually orthogonal system.
Definition: Vec.cpp:41
CylMesh::setY0
void setY0(const Eref &e, double v)
Definition: CylMesh.cpp:251
CylMesh::setZ1
void setZ1(const Eref &e, double v)
Definition: CylMesh.cpp:317
CylMesh::getDiffLength
double getDiffLength(const Eref &e) const
Definition: CylMesh.cpp:400
Field::get
static A get(const ObjId &dest, const string &field)
Definition: SetGet.h:284
CubeMesh::setPreserveNumEntries
void setPreserveNumEntries(bool v)
Definition: CubeMesh.cpp:581
testCubeMeshMultiJunctionTwoD
void testCubeMeshMultiJunctionTwoD()
Definition: testMesh.cpp:1636
CubeMesh::ABUTX
static const unsigned int ABUTX
Definition: CubeMesh.h:287
CubeMesh::getMeshDimensions
unsigned int getMeshDimensions(unsigned int fid) const
Virtual function to return dimensions of specified entry.
Definition: CubeMesh.cpp:778
CubeMesh::numDims
unsigned int numDims() const
Utility function for returning # of dimensions in mesh.
Definition: CubeMesh.cpp:1148
CubeMesh::setSurface
void setSurface(vector< unsigned int > v)
Definition: CubeMesh.cpp:677
DataId
unsigned int DataId
Definition: header.h:47
testCubeMeshExtendStencil
void testCubeMeshExtendStencil()
Definition: testMesh.cpp:726
CubeMesh::getDy
double getDy() const
Definition: CubeMesh.cpp:519
CubeMesh::setX0
void setX0(double v)
Definition: CubeMesh.cpp:435
CylMesh::getZ1
double getZ1(const Eref &e) const
Definition: CylMesh.cpp:325
CubeMesh::innerGetNumEntries
unsigned int innerGetNumEntries() const
Definition: CubeMesh.cpp:891
testCubeMesh
void testCubeMesh()
Definition: testMesh.cpp:600
Shell
Definition: Shell.h:43
CylMesh::selectGridVolume
double selectGridVolume(double h) const
Definition: CylMesh.cpp:940
CylMesh::innerSetCoords
void innerSetCoords(const Eref &e, const vector< double > &v)
Definition: CylMesh.cpp:343
SetGet2::set
static bool set(const ObjId &dest, const string &field, A1 arg1, A2 arg2)
Definition: SetGet.h:365