Skip to content

simplified dc functions #369

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: development
Choose a base branch
from
Open

simplified dc functions #369

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9cb3fd3
simplified dc functions
drbergman Apr 1, 2025
d3d5d92
check size of values in (un)fix dcs
drbergman Apr 1, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
249 changes: 206 additions & 43 deletions BioFVM/BioFVM_microenvironment.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@
#include "BioFVM_solvers.h"
#include "BioFVM_vector.h"
#include <cmath>
#include <algorithm>

#include "BioFVM_basic_agent.h"

Expand Down Expand Up @@ -185,15 +186,175 @@ Microenvironment::Microenvironment(std::string name)
return;
}

void Microenvironment::fix_substrate_at_voxel( std::string substrate, int voxel_index, double new_value )
{
int substrate_index = find_existing_density_index( substrate );
return fix_substrate_at_voxel(substrate_index, voxel_index, new_value);
}

void Microenvironment::fix_substrate_at_voxel( std::string substrate, int voxel_index )
{
int substrate_index = find_existing_density_index( substrate );
return fix_substrate_at_voxel(substrate_index, voxel_index);
}

void Microenvironment::fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices, double new_value )
{
int substrate_index = find_existing_density_index( substrate );
return fix_substrate_at_voxels(substrate_index, voxel_indices, new_value);
}

void Microenvironment::fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices, std::vector<double>& new_values )
{
int substrate_index = find_existing_density_index( substrate );
return fix_substrate_at_voxels(substrate_index, voxel_indices, new_values);
}

void Microenvironment::fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices )
{
int substrate_index = find_existing_density_index( substrate );
return fix_substrate_at_voxels(substrate_index, voxel_indices);
}

void Microenvironment::fix_substrate_at_voxel( int substrate_index, int voxel_index, double new_value )
{
dirichlet_value_vectors[voxel_index][substrate_index] = new_value;
fix_substrate_at_voxel(substrate_index, voxel_index);
}

void Microenvironment::fix_substrate_at_voxel( int substrate_index, int voxel_index )
{
dirichlet_activation_vectors[voxel_index][substrate_index] = true;
mesh.voxels[voxel_index].is_Dirichlet = true;
dirichlet_activation_vector[substrate_index] = true;
return;
}

void Microenvironment::fix_substrates_at_voxel( int voxel_index , std::vector<double>& new_values )
{
if (new_values.size() != dirichlet_value_vectors[voxel_index].size())
{
std::cerr << "Error: Incorrect number of values passed in to fix_substrates_at_voxel. Expected " << dirichlet_value_vectors[voxel_index].size() << " values, but got " << new_values.size() << "." << std::endl;
return;
}
dirichlet_value_vectors[voxel_index] = new_values;
return fix_substrates_at_voxel( voxel_index );
}

void Microenvironment::fix_substrates_at_voxel( int voxel_index )
{
for( unsigned int i=0 ; i < dirichlet_activation_vectors[voxel_index].size() ; i++ )
{
fix_substrate_at_voxel( i , voxel_index );
}
return;
}

void Microenvironment::fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices, double new_value )
{
for( auto voxel_index : voxel_indices )
{
dirichlet_value_vectors[voxel_index][substrate_index] = new_value;
}
return fix_substrate_at_voxels( substrate_index, voxel_indices );
}

void Microenvironment::fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices, std::vector<double>& new_values )
{
if (new_values.size() != voxel_indices.size())
{
std::cerr << "Error: new_values size (" << new_values.size() << ") does not match voxel_indices size (" << voxel_indices.size() << ") in Microenvironment::fix_substrate_at_voxels" << std::endl;
return;
}
for( unsigned int i=0 ; i < voxel_indices.size() ; i++ )
{
dirichlet_value_vectors[voxel_indices[i]][substrate_index] = new_values[i];
}
return fix_substrate_at_voxels( substrate_index, voxel_indices );
}

void Microenvironment::fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices )
{
for( auto voxel_index : voxel_indices )
{
dirichlet_activation_vectors[voxel_index][substrate_index] = true;
mesh.voxels[voxel_index].is_Dirichlet = true;
}

// do not allow dirichlet_activation_vector[substrate_index] to be set to true if no voxels are changed
if (voxel_indices.size() > 0)
{
dirichlet_activation_vector[substrate_index] = true;
}
return;
}

void Microenvironment::unfix_substrate_at_voxel( std::string substrate, int voxel_index )
{
int substrate_index = find_existing_density_index( substrate );
return unfix_substrate_at_voxel(substrate_index, voxel_index);
}

void Microenvironment::unfix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices )
{
int substrate_index = find_existing_density_index( substrate );
return unfix_substrate_at_voxels(substrate_index, voxel_indices);
}

void Microenvironment::unfix_substrate_at_voxel( int substrate_index, int voxel_index )
{
dirichlet_activation_vectors[voxel_index][substrate_index] = false;
if (std::find(dirichlet_activation_vectors[voxel_index].begin(), dirichlet_activation_vectors[voxel_index].end(), true) == dirichlet_activation_vectors[voxel_index].end())
{
mesh.voxels[voxel_index].is_Dirichlet = false;
}
sync_substrate_dirichlet_activation(substrate_index);
return;
}

void Microenvironment::unfix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices )
{
// 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
for( auto voxel_index : voxel_indices )
{
dirichlet_activation_vectors[voxel_index][substrate_index] = false;
if (std::find(dirichlet_activation_vectors[voxel_index].begin(), dirichlet_activation_vectors[voxel_index].end(), true) == dirichlet_activation_vectors[voxel_index].end())
{
mesh.voxels[voxel_index].is_Dirichlet = false;
}
}
sync_substrate_dirichlet_activation(substrate_index);
}

void Microenvironment::unfix_substrates_at_voxel( int voxel_index )
{
for( unsigned int i=0 ; i < dirichlet_activation_vectors[voxel_index].size() ; i++ )
{
dirichlet_activation_vectors[voxel_index][i] = false;
sync_substrate_dirichlet_activation(i);
}
mesh.voxels[voxel_index].is_Dirichlet = false;
return;
}

void Microenvironment::sync_substrate_dirichlet_activation( int substrate_index )
{
for (auto voxel_activation_vector : dirichlet_activation_vectors)
{
if (voxel_activation_vector[substrate_index])
{
dirichlet_activation_vector[substrate_index] = true;
return;
}
}
dirichlet_activation_vector[substrate_index] = false;
return;
}

void Microenvironment::add_dirichlet_node( int voxel_index, std::vector<double>& value )
{
mesh.voxels[voxel_index].is_Dirichlet=true;
/*
dirichlet_indices.push_back( voxel_index );
dirichlet_value_vectors.push_back( value );
*/

dirichlet_value_vectors[voxel_index] = value; // .assign( mesh.voxels.size(), one );
dirichlet_value_vectors[voxel_index] = value;

return;
}
Expand Down Expand Up @@ -569,6 +730,15 @@ int Microenvironment::find_density_index( std::string name )
return -1;
}

int Microenvironment::find_existing_density_index( std::string name )
{
int i = find_density_index( name );
if( i != -1 )
{ return i; }
std::cerr << "Error: density named " << name << " not found." << std::endl;
exit(-1);
}

void Microenvironment::set_density( int index , std::string name , std::string units )
{
// fix in PhysiCell preview November 2017
Expand Down Expand Up @@ -1282,13 +1452,11 @@ void set_microenvironment_initial_condition( void )
// set Dirichlet conditions along the xmin outer edges
for( unsigned int j=0 ; j < microenvironment.mesh.y_coordinates.size() ; j++ )
{
// set the value
microenvironment.add_dirichlet_node( microenvironment.voxel_index(I,j,k) , default_microenvironment_options.Dirichlet_xmin_values );

// set the activation
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(I,j,k) ,
default_microenvironment_options.Dirichlet_xmin );

for ( size_t n = 0; n < default_microenvironment_options.Dirichlet_xmin_values.size(); n++)
{
if (!default_microenvironment_options.Dirichlet_xmin[n]) { continue; }
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(I,j,k), default_microenvironment_options.Dirichlet_xmin_values[n] );
}
}
}
}
Expand All @@ -1302,12 +1470,11 @@ void set_microenvironment_initial_condition( void )
// set Dirichlet conditions along the xmax outer edges
for( unsigned int j=0 ; j < microenvironment.mesh.y_coordinates.size() ; j++ )
{
// set the values
microenvironment.add_dirichlet_node( microenvironment.voxel_index(I,j,k) , default_microenvironment_options.Dirichlet_xmax_values );

// set the activation
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(I,j,k) ,
default_microenvironment_options.Dirichlet_xmax );
for (size_t n = 0; n < default_microenvironment_options.Dirichlet_xmax_values.size(); n++)
{
if (!default_microenvironment_options.Dirichlet_xmax[n]) { continue; }
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(I,j,k), default_microenvironment_options.Dirichlet_xmax_values[n] );
}
}
}
}
Expand All @@ -1321,12 +1488,11 @@ void set_microenvironment_initial_condition( void )
// set Dirichlet conditions along the ymin outer edges
for( unsigned int i=0 ; i < microenvironment.mesh.x_coordinates.size() ; i++ )
{
// set the values
microenvironment.add_dirichlet_node( microenvironment.voxel_index(i,J,k) , default_microenvironment_options.Dirichlet_ymin_values );

// set the activation
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(i,J,k) ,
default_microenvironment_options.Dirichlet_ymin );
for (size_t n = 0; n < default_microenvironment_options.Dirichlet_ymin_values.size(); n++)
{
if (!default_microenvironment_options.Dirichlet_ymin[n]) { continue; }
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(i,J,k), default_microenvironment_options.Dirichlet_ymin_values[n] );
}
}
}
}
Expand All @@ -1340,12 +1506,11 @@ void set_microenvironment_initial_condition( void )
// set Dirichlet conditions along the ymin outer edges
for( unsigned int i=0 ; i < microenvironment.mesh.x_coordinates.size() ; i++ )
{
// set the value
microenvironment.add_dirichlet_node( microenvironment.voxel_index(i,J,k) , default_microenvironment_options.Dirichlet_ymax_values );

// set the activation
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(i,J,k) ,
default_microenvironment_options.Dirichlet_ymax );
for (size_t n = 0; n < default_microenvironment_options.Dirichlet_ymax_values.size(); n++)
{
if (!default_microenvironment_options.Dirichlet_ymax[n]) { continue; }
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(i,J,k), default_microenvironment_options.Dirichlet_ymax_values[n] );
}
}
}
}
Expand All @@ -1362,12 +1527,11 @@ void set_microenvironment_initial_condition( void )
// set Dirichlet conditions along the ymin outer edges
for( unsigned int i=0 ; i < microenvironment.mesh.x_coordinates.size() ; i++ )
{
// set the value
microenvironment.add_dirichlet_node( microenvironment.voxel_index(i,j,K) , default_microenvironment_options.Dirichlet_zmin_values );

// set the activation
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(i,j,K) ,
default_microenvironment_options.Dirichlet_zmin );
for ( size_t n = 0; n < default_microenvironment_options.Dirichlet_zmin_values.size(); n++)
{
if (!default_microenvironment_options.Dirichlet_zmin[n]) { continue; }
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(i,j,K), default_microenvironment_options.Dirichlet_zmin_values[n] );
}
}
}
}
Expand All @@ -1381,12 +1545,11 @@ void set_microenvironment_initial_condition( void )
// set Dirichlet conditions along the ymin outer edges
for( unsigned int i=0 ; i < microenvironment.mesh.x_coordinates.size() ; i++ )
{
// set the value
microenvironment.add_dirichlet_node( microenvironment.voxel_index(i,j,K) , default_microenvironment_options.Dirichlet_zmax_values );

// set the activation
microenvironment.set_substrate_dirichlet_activation( microenvironment.voxel_index(i,j,K) ,
default_microenvironment_options.Dirichlet_zmax );
for ( size_t n = 0; n < default_microenvironment_options.Dirichlet_zmax_values.size(); n++)
{
if (!default_microenvironment_options.Dirichlet_zmax[n]) { continue; }
microenvironment.fix_substrate_at_voxel( n, microenvironment.voxel_index(i,j,K), default_microenvironment_options.Dirichlet_zmax_values[n] );
}
}
}
}
Expand Down
34 changes: 24 additions & 10 deletions BioFVM/BioFVM_microenvironment.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,18 +121,10 @@ class Microenvironment

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

/*
std::vector<int> dirichlet_indices;
std::vector< std::vector<double> > dirichlet_value_vectors;
std::vector<bool> dirichlet_node_map;
*/
std::vector< std::vector<double> > dirichlet_value_vectors;
std::vector<bool> dirichlet_activation_vector;

/* new in Version 1.7.0 -- activation vectors can be specified
on a voxel-by-voxel basis */

std::vector< std::vector<bool> > dirichlet_activation_vectors;
std::vector<bool> dirichlet_activation_vector; // whether a given substrate has a DC set somewhere

public:

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

int find_density_index( std::string name );
int find_existing_density_index( std::string name );

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

void display_information( std::ostream& os );


void fix_substrate_at_voxel( std::string substrate, int voxel_index, double new_value );
void fix_substrate_at_voxel( std::string substrate, int voxel_index );
void fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices, double new_value );
void fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices, std::vector<double>& new_values );
void fix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices );
void fix_substrate_at_voxel( int substrate_index, int voxel_index, double new_value );
void fix_substrate_at_voxel( int substrate_index, int voxel_index );
void fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices, double new_value );
void fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices, std::vector<double>& new_values );
void fix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices );
void fix_substrates_at_voxel( int voxel_index, std::vector<double>& new_values );
void fix_substrates_at_voxel( int voxel_index );

void unfix_substrate_at_voxel( std::string substrate, int voxel_index );
void unfix_substrate_at_voxels( std::string substrate, std::vector<int>& voxel_indices );
void unfix_substrate_at_voxel( int substrate_index, int voxel_index );
void unfix_substrate_at_voxels( int substrate_index, std::vector<int>& voxel_indices );
void unfix_substrates_at_voxel( int voxel_index );

void sync_substrate_dirichlet_activation( int substrate_index );

void add_dirichlet_node( int voxel_index, std::vector<double>& value );
void update_dirichlet_node( int voxel_index , std::vector<double>& new_value );
void update_dirichlet_node( int voxel_index , int substrate_index , double new_value );
Expand Down
Loading