Skip to content

Commit 9cb3fd3

Browse files
drbergmandrbergman
committed
simplified dc functions
1 parent cbf625e commit 9cb3fd3

File tree

2 files changed

+225
-53
lines changed

2 files changed

+225
-53
lines changed

‎BioFVM/BioFVM_microenvironment.cpp‎

Lines changed: 201 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,7 @@
5050
#include "BioFVM_solvers.h"
5151
#include "BioFVM_vector.h"
5252
#include <cmath>
53+
#include <algorithm>
5354

5455
#include "BioFVM_basic_agent.h"
5556

@@ -185,15 +186,170 @@ Microenvironment::Microenvironment(std::string name)
185186
return;
186187
}
187188

189+
void Microenvironment::fix_substrate_at_voxel( std::string substrate, int voxel_index, double new_value )
190+
{
191+
int substrate_index = find_existing_density_index( substrate );
192+
return fix_substrate_at_voxel(substrate_index, voxel_index, new_value);
193+
}
194+
195+
void Microenvironment::fix_substrate_at_voxel( std::string substrate, int voxel_index )
196+
{
197+
int substrate_index = find_existing_density_index( substrate );
198+
return fix_substrate_at_voxel(substrate_index, voxel_index);
199+
}
200+
201+
void Microenvironment::fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices, double new_value )
202+
{
203+
int substrate_index = find_existing_density_index( substrate );
204+
return fix_substrate_at_voxels(substrate_index, voxel_indices, new_value);
205+
}
206+
207+
void Microenvironment::fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices, std::vector<double>& new_values )
208+
{
209+
int substrate_index = find_existing_density_index( substrate );
210+
return fix_substrate_at_voxels(substrate_index, voxel_indices, new_values);
211+
}
212+
213+
void Microenvironment::fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices )
214+
{
215+
int substrate_index = find_existing_density_index( substrate );
216+
return fix_substrate_at_voxels(substrate_index, voxel_indices);
217+
}
218+
219+
void Microenvironment::fix_substrate_at_voxel( int substrate_index, int voxel_index, double new_value )
220+
{
221+
dirichlet_value_vectors[voxel_index][substrate_index] = new_value;
222+
fix_substrate_at_voxel(substrate_index, voxel_index);
223+
}
224+
225+
void Microenvironment::fix_substrate_at_voxel( int substrate_index, int voxel_index )
226+
{
227+
dirichlet_activation_vectors[voxel_index][substrate_index] = true;
228+
mesh.voxels[voxel_index].is_Dirichlet = true;
229+
dirichlet_activation_vector[substrate_index] = true;
230+
return;
231+
}
232+
233+
void Microenvironment::fix_substrates_at_voxel( int voxel_index , std::vector<double>& new_values )
234+
{
235+
dirichlet_value_vectors[voxel_index] = new_values;
236+
return fix_substrates_at_voxel( voxel_index );
237+
}
238+
239+
void Microenvironment::fix_substrates_at_voxel( int voxel_index )
240+
{
241+
for( unsigned int i=0 ; i < dirichlet_activation_vectors[voxel_index].size() ; i++ )
242+
{
243+
fix_substrate_at_voxel( i , voxel_index );
244+
}
245+
return;
246+
}
247+
248+
void Microenvironment::fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices, double new_value )
249+
{
250+
for( auto voxel_index : voxel_indices )
251+
{
252+
dirichlet_value_vectors[voxel_index][substrate_index] = new_value;
253+
}
254+
return fix_substrate_at_voxels( substrate_index, voxel_indices );
255+
}
256+
257+
void Microenvironment::fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices, std::vector<double>& new_values )
258+
{
259+
if (new_values.size() != voxel_indices.size())
260+
{
261+
std::cerr << "Error: new_values size does not match voxel_indices size in Microenvironment::fix_substrate_at_voxels" << std::endl;
262+
return;
263+
}
264+
for( unsigned int i=0 ; i < voxel_indices.size() ; i++ )
265+
{
266+
dirichlet_value_vectors[voxel_indices[i]][substrate_index] = new_values[i];
267+
}
268+
return fix_substrate_at_voxels( substrate_index, voxel_indices );
269+
}
270+
271+
void Microenvironment::fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices )
272+
{
273+
for( auto voxel_index : voxel_indices )
274+
{
275+
dirichlet_activation_vectors[voxel_index][substrate_index] = true;
276+
mesh.voxels[voxel_index].is_Dirichlet = true;
277+
}
278+
279+
// do not allow dirichlet_activation_vector[substrate_index] to be set to true if no voxels are changed
280+
if (voxel_indices.size() > 0)
281+
{
282+
dirichlet_activation_vector[substrate_index] = true;
283+
}
284+
return;
285+
}
286+
287+
void Microenvironment::unfix_substrate_at_voxel( std::string substrate, int voxel_index )
288+
{
289+
int substrate_index = find_existing_density_index( substrate );
290+
return unfix_substrate_at_voxel(substrate_index, voxel_index);
291+
}
292+
293+
void Microenvironment::unfix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices )
294+
{
295+
int substrate_index = find_existing_density_index( substrate );
296+
return unfix_substrate_at_voxels(substrate_index, voxel_indices);
297+
}
298+
299+
void Microenvironment::unfix_substrate_at_voxel( int substrate_index, int voxel_index )
300+
{
301+
dirichlet_activation_vectors[voxel_index][substrate_index] = false;
302+
if (std::find(dirichlet_activation_vectors[voxel_index].begin(), dirichlet_activation_vectors[voxel_index].end(), true) == dirichlet_activation_vectors[voxel_index].end())
303+
{
304+
mesh.voxels[voxel_index].is_Dirichlet = false;
305+
}
306+
sync_substrate_dirichlet_activation(substrate_index);
307+
return;
308+
}
309+
310+
void Microenvironment::unfix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices )
311+
{
312+
// this allows for unfixing a substrate at multiple voxels at once without checking if the substrate remains a DC substrate after removing from each voxel
313+
for( auto voxel_index : voxel_indices )
314+
{
315+
dirichlet_activation_vectors[voxel_index][substrate_index] = false;
316+
if (std::find(dirichlet_activation_vectors[voxel_index].begin(), dirichlet_activation_vectors[voxel_index].end(), true) == dirichlet_activation_vectors[voxel_index].end())
317+
{
318+
mesh.voxels[voxel_index].is_Dirichlet = false;
319+
}
320+
}
321+
sync_substrate_dirichlet_activation(substrate_index);
322+
}
323+
324+
void Microenvironment::unfix_substrates_at_voxel( int voxel_index )
325+
{
326+
for( unsigned int i=0 ; i < dirichlet_activation_vectors[voxel_index].size() ; i++ )
327+
{
328+
dirichlet_activation_vectors[voxel_index][i] = false;
329+
sync_substrate_dirichlet_activation(i);
330+
}
331+
mesh.voxels[voxel_index].is_Dirichlet = false;
332+
return;
333+
}
334+
335+
void Microenvironment::sync_substrate_dirichlet_activation( int substrate_index )
336+
{
337+
for (auto voxel_activation_vector : dirichlet_activation_vectors)
338+
{
339+
if (voxel_activation_vector[substrate_index])
340+
{
341+
dirichlet_activation_vector[substrate_index] = true;
342+
return;
343+
}
344+
}
345+
dirichlet_activation_vector[substrate_index] = false;
346+
return;
347+
}
348+
188349
void Microenvironment::add_dirichlet_node( int voxel_index, std::vector<double>& value )
189350
{
190351
mesh.voxels[voxel_index].is_Dirichlet=true;
191-
/*
192-
dirichlet_indices.push_back( voxel_index );
193-
dirichlet_value_vectors.push_back( value );
194-
*/
195-
196-
dirichlet_value_vectors[voxel_index] = value; // .assign( mesh.voxels.size(), one );
352+
dirichlet_value_vectors[voxel_index] = value;
197353

198354
return;
199355
}
@@ -569,6 +725,15 @@ int Microenvironment::find_density_index( std::string name )
569725
return -1;
570726
}
571727

728+
int Microenvironment::find_existing_density_index( std::string name )
729+
{
730+
int i = find_density_index( name );
731+
if( i != -1 )
732+
{ return i; }
733+
std::cout << "Error: density named " << name << " not found." << std::endl;
734+
exit(-1);
735+
}
736+
572737
void Microenvironment::set_density( int index , std::string name , std::string units )
573738
{
574739
// fix in PhysiCell preview November 2017
@@ -1282,13 +1447,11 @@ void set_microenvironment_initial_condition( void )
12821447
// set Dirichlet conditions along the xmin outer edges
12831448
for( unsigned int j=0 ; j < microenvironment.mesh.y_coordinates.size() ; j++ )
12841449
{
1285-
// set the value
1286-
microenvironment.add_dirichlet_node( microenvironment.voxel_index(I,j,k) , default_microenvironment_options.Dirichlet_xmin_values );
1287-
1288-
// set the activation
1289-
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(I,j,k) ,
1290-
default_microenvironment_options.Dirichlet_xmin );
1291-
1450+
for ( size_t n = 0; n < default_microenvironment_options.Dirichlet_xmin_values.size(); n++)
1451+
{
1452+
if (!default_microenvironment_options.Dirichlet_xmin[n]) { continue; }
1453+
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(I,j,k), default_microenvironment_options.Dirichlet_xmin_values[n] );
1454+
}
12921455
}
12931456
}
12941457
}
@@ -1302,12 +1465,11 @@ void set_microenvironment_initial_condition( void )
13021465
// set Dirichlet conditions along the xmax outer edges
13031466
for( unsigned int j=0 ; j < microenvironment.mesh.y_coordinates.size() ; j++ )
13041467
{
1305-
// set the values
1306-
microenvironment.add_dirichlet_node( microenvironment.voxel_index(I,j,k) , default_microenvironment_options.Dirichlet_xmax_values );
1307-
1308-
// set the activation
1309-
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(I,j,k) ,
1310-
default_microenvironment_options.Dirichlet_xmax );
1468+
for (size_t n = 0; n < default_microenvironment_options.Dirichlet_xmax_values.size(); n++)
1469+
{
1470+
if (!default_microenvironment_options.Dirichlet_xmax[n]) { continue; }
1471+
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(I,j,k), default_microenvironment_options.Dirichlet_xmax_values[n] );
1472+
}
13111473
}
13121474
}
13131475
}
@@ -1321,12 +1483,11 @@ void set_microenvironment_initial_condition( void )
13211483
// set Dirichlet conditions along the ymin outer edges
13221484
for( unsigned int i=0 ; i < microenvironment.mesh.x_coordinates.size() ; i++ )
13231485
{
1324-
// set the values
1325-
microenvironment.add_dirichlet_node( microenvironment.voxel_index(i,J,k) , default_microenvironment_options.Dirichlet_ymin_values );
1326-
1327-
// set the activation
1328-
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(i,J,k) ,
1329-
default_microenvironment_options.Dirichlet_ymin );
1486+
for (size_t n = 0; n < default_microenvironment_options.Dirichlet_ymin_values.size(); n++)
1487+
{
1488+
if (!default_microenvironment_options.Dirichlet_ymin[n]) { continue; }
1489+
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(i,J,k), default_microenvironment_options.Dirichlet_ymin_values[n] );
1490+
}
13301491
}
13311492
}
13321493
}
@@ -1340,12 +1501,11 @@ void set_microenvironment_initial_condition( void )
13401501
// set Dirichlet conditions along the ymin outer edges
13411502
for( unsigned int i=0 ; i < microenvironment.mesh.x_coordinates.size() ; i++ )
13421503
{
1343-
// set the value
1344-
microenvironment.add_dirichlet_node( microenvironment.voxel_index(i,J,k) , default_microenvironment_options.Dirichlet_ymax_values );
1345-
1346-
// set the activation
1347-
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(i,J,k) ,
1348-
default_microenvironment_options.Dirichlet_ymax );
1504+
for (size_t n = 0; n < default_microenvironment_options.Dirichlet_ymax_values.size(); n++)
1505+
{
1506+
if (!default_microenvironment_options.Dirichlet_ymax[n]) { continue; }
1507+
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(i,J,k), default_microenvironment_options.Dirichlet_ymax_values[n] );
1508+
}
13491509
}
13501510
}
13511511
}
@@ -1362,12 +1522,11 @@ void set_microenvironment_initial_condition( void )
13621522
// set Dirichlet conditions along the ymin outer edges
13631523
for( unsigned int i=0 ; i < microenvironment.mesh.x_coordinates.size() ; i++ )
13641524
{
1365-
// set the value
1366-
microenvironment.add_dirichlet_node( microenvironment.voxel_index(i,j,K) , default_microenvironment_options.Dirichlet_zmin_values );
1367-
1368-
// set the activation
1369-
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(i,j,K) ,
1370-
default_microenvironment_options.Dirichlet_zmin );
1525+
for ( size_t n = 0; n < default_microenvironment_options.Dirichlet_zmin_values.size(); n++)
1526+
{
1527+
if (!default_microenvironment_options.Dirichlet_zmin[n]) { continue; }
1528+
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(i,j,K), default_microenvironment_options.Dirichlet_zmin_values[n] );
1529+
}
13711530
}
13721531
}
13731532
}
@@ -1381,12 +1540,11 @@ void set_microenvironment_initial_condition( void )
13811540
// set Dirichlet conditions along the ymin outer edges
13821541
for( unsigned int i=0 ; i < microenvironment.mesh.x_coordinates.size() ; i++ )
13831542
{
1384-
// set the value
1385-
microenvironment.add_dirichlet_node( microenvironment.voxel_index(i,j,K) , default_microenvironment_options.Dirichlet_zmax_values );
1386-
1387-
// set the activation
1388-
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(i,j,K) ,
1389-
default_microenvironment_options.Dirichlet_zmax );
1543+
for ( size_t n = 0; n < default_microenvironment_options.Dirichlet_zmax_values.size(); n++)
1544+
{
1545+
if (!default_microenvironment_options.Dirichlet_zmax[n]) { continue; }
1546+
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(i,j,K), default_microenvironment_options.Dirichlet_zmax_values[n] );
1547+
}
13901548
}
13911549
}
13921550
}

‎BioFVM/BioFVM_microenvironment.h‎

Lines changed: 24 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -121,18 +121,10 @@ class Microenvironment
121121

122122
// on "resize density" type operations, need to extend all of these
123123

124-
/*
125-
std::vector<int> dirichlet_indices;
126124
std::vector< std::vector<double> > dirichlet_value_vectors;
127-
std::vector<bool> dirichlet_node_map;
128-
*/
129-
std::vector< std::vector<double> > dirichlet_value_vectors;
130-
std::vector<bool> dirichlet_activation_vector;
131125

132-
/* new in Version 1.7.0 -- activation vectors can be specified
133-
on a voxel-by-voxel basis */
134-
135126
std::vector< std::vector<bool> > dirichlet_activation_vectors;
127+
std::vector<bool> dirichlet_activation_vector; // whether a given substrate has a DC set somewhere
136128

137129
public:
138130

@@ -190,6 +182,7 @@ class Microenvironment
190182
void set_density( int index , std::string name , std::string units , double diffusion_constant , double decay_rate );
191183

192184
int find_density_index( std::string name );
185+
int find_existing_density_index( std::string name );
193186

194187
int voxel_index( int i, int j, int k );
195188
std::vector<unsigned int> cartesian_indices( int n );
@@ -237,7 +230,28 @@ class Microenvironment
237230
void simulate_cell_sources_and_sinks( double dt );
238231

239232
void display_information( std::ostream& os );
240-
233+
234+
void fix_substrate_at_voxel( std::string substrate, int voxel_index, double new_value );
235+
void fix_substrate_at_voxel( std::string substrate, int voxel_index );
236+
void fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices, double new_value );
237+
void fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices, std::vector<double>& new_values );
238+
void fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices );
239+
void fix_substrate_at_voxel( int substrate_index, int voxel_index, double new_value );
240+
void fix_substrate_at_voxel( int substrate_index, int voxel_index );
241+
void fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices, double new_value );
242+
void fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices, std::vector<double>& new_values );
243+
void fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices );
244+
void fix_substrates_at_voxel( int voxel_ind, std::vector<double>& new_values );
245+
void fix_substrates_at_voxel( int voxel_ind );
246+
247+
void unfix_substrate_at_voxel( std::string substrate, int voxel_index );
248+
void unfix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices );
249+
void unfix_substrate_at_voxel( int substrate_index, int voxel_index );
250+
void unfix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices );
251+
void unfix_substrates_at_voxel( int voxel_index );
252+
253+
void sync_substrate_dirichlet_activation( int substrate_index );
254+
241255
void add_dirichlet_node( int voxel_index, std::vector<double>& value );
242256
void update_dirichlet_node( int voxel_index , std::vector<double>& new_value );
243257
void update_dirichlet_node( int voxel_index , int substrate_index , double new_value );

0 commit comments

Comments
 (0)