46 std::vector<double>& rStimulusMagnitudes,
47 std::vector<double>& rDurations,
48 std::vector<double>& rPeriods,
49 std::vector<double>& rStarts,
50 std::vector<double>& rEnds)
51 : mrElectrodePairs(rElectrodePairs),
52 mrMagnitudes(rStimulusMagnitudes),
53 mrDurations(rDurations),
57 mGroundSecondElectrode(false)
59 if ((rElectrodePairs.size() != rStimulusMagnitudes.size()) ||
60 (rElectrodePairs.size() != rDurations.size()) ||
61 (rElectrodePairs.size() != rPeriods.size()) ||
62 (rElectrodePairs.size() != rStarts.size()) ||
63 (rElectrodePairs.size() != rEnds.size()))
65 EXCEPTION (
"Vector of electrode pairs and vector of stimulation parameters must have the same size");
80 std::vector<unsigned> nodes_in_all_electrodes;
81 for (
unsigned global_node_index = 0; global_node_index < this->mpMesh->GetNumNodes(); global_node_index++)
83 if (this->mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(global_node_index))
85 for (
unsigned pair_index = 0; pair_index <mrElectrodePairs.size(); pair_index++)
87 if (mrElectrodePairs[pair_index].first->DoesContain( this->mpMesh->GetNode(global_node_index)->GetPoint()))
89 nodes_in_all_electrodes.push_back( this->mpMesh->GetNode(global_node_index)->GetIndex() );
91 if (mrElectrodePairs[pair_index].second->DoesContain( this->mpMesh->GetNode(global_node_index)->GetPoint()))
93 nodes_in_all_electrodes.push_back( this->mpMesh->GetNode(global_node_index)->GetIndex() );
99 for (
unsigned node_index = 0; node_index < nodes_in_all_electrodes.size(); node_index++)
101 unsigned number_of_hits = 0;
102 for (
unsigned node_to_check = 0; node_to_check < nodes_in_all_electrodes.size(); node_to_check++)
104 if (nodes_in_all_electrodes[node_index] == nodes_in_all_electrodes[node_to_check] )
109 if (number_of_hits>1)
111 EXCEPTION(
"Two or more electrodes intersect with each other");
125 assert(this->mpMesh!=NULL);
128 CheckForElectrodesIntersection();
137 for (
unsigned pair_index = 0; pair_index < mrElectrodePairs.size(); pair_index++)
140 if (!mGroundSecondElectrode)
143 double flux_electrode_1 = ComputeElectrodeTotalFlux(mrElectrodePairs[pair_index].first, mMagnitudesElectrode1[pair_index]);
144 double flux_electrode_2 = ComputeElectrodeTotalFlux(mrElectrodePairs[pair_index].second, mMagnitudesElectrode2[pair_index]);
147 mMagnitudesElectrode2[pair_index] = -mMagnitudesElectrode1[pair_index]*flux_electrode_1/flux_electrode_2;
150 assert( flux_electrode_2 != 0.0);
151 assert( ! std::isnan(mMagnitudesElectrode2[pair_index]));
156 this->mGroundedRegions.push_back( mrElectrodePairs[pair_index].second );
164 boost::shared_ptr<RegularStimulus> p_stimulus;
165 for (
unsigned pair_index = 0; pair_index < mrElectrodePairs.size(); pair_index++)
167 if (mrElectrodePairs[pair_index].first->DoesContain(pNode->
GetPoint()) )
169 p_stimulus.reset(
new RegularStimulus(mMagnitudesElectrode1[pair_index], mrDurations[pair_index], mrPeriods[pair_index], mrStarts[pair_index], mrEnds[pair_index]));
172 else if (mrElectrodePairs[pair_index].second->DoesContain(pNode->
GetPoint()) )
174 p_stimulus.reset (
new RegularStimulus(mMagnitudesElectrode2[pair_index], mrDurations[pair_index], mrPeriods[pair_index], mrStarts[pair_index], mrEnds[pair_index]));
179 p_stimulus.reset (
new RegularStimulus(0.0, mrDurations[pair_index], mrPeriods[pair_index], mrStarts[pair_index], mrEnds[pair_index]));
191 c_vector<double, DIM+1> phi;
193 double total_electrode_flux = 0.0;
197 node_iter != this->mpMesh->GetNodeIteratorEnd();
202 unsigned node_index = node_iter->GetIndex();
203 assert(node_index < this->mpMesh->GetNumNodes());
204 double contribution_of_this_node = 0.0;
206 for (std::set<unsigned>::iterator iter = this->mpMesh->GetNode(node_index)->rGetContainingElementIndices().begin();
207 iter != this->mpMesh->GetNode(node_index)->rGetContainingElementIndices().end();
213 c_matrix<double, DIM, DIM> jacobian;
214 c_matrix<double, DIM, DIM> inverse_jacobian;
215 double jacobian_determinant;
216 this->mpMesh->GetInverseJacobianForElement(p_element->
GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
218 double contribution_of_this_element = 0.0;
220 for (
unsigned quad_index=0; quad_index < pQuadRule->
GetNumQuadPoints(); quad_index++)
223 BasisFunction::ComputeBasisFunctions(quad_point, phi);
225 double interpolated_stimulus = 0.0;
227 for (
unsigned node_index_in_element = 0; node_index_in_element < p_element->
GetNumNodes(); node_index_in_element++)
231 interpolated_stimulus += stimulusMagnitude*phi(node_index_in_element);
232 contribution_of_this_element += interpolated_stimulus * phi(node_index_in_element) * jacobian_determinant * pQuadRule->
GetWeight(quad_index);
236 contribution_of_this_node += contribution_of_this_element;
239 total_electrode_flux += contribution_of_this_node;
243 int mpi_ret = MPI_Allreduce(&total_electrode_flux, &ret, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
244 assert(mpi_ret == MPI_SUCCESS);
246 MPI_Allreduce(&total_electrode_flux, &ret, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
252 assert(ret < DBL_MAX);
ElectrodesStimulusFactory(std::vector< std::pair< AbstractChasteRegion< DIM > *, AbstractChasteRegion< DIM > * > > &rElectrodePairs, std::vector< double > &rStimulusMagnitudes, std::vector< double > &rDurations, std::vector< double > &rPeriods, std::vector< double > &rStarts, std::vector< double > &rEnds)