42 assert(this->mpLinearSystem == NULL);
50 std::cout <<
"Using absolute tolerance: " << mpConfig->GetAbsoluteTolerance() <<
"\n";
52 this->mpLinearSystem->SetAbsoluteTolerance(mpConfig->GetAbsoluteTolerance());
57 std::cout <<
"Using relative tolerance: " << mpConfig->GetRelativeTolerance() <<
"\n";
59 this->mpLinearSystem->SetRelativeTolerance(mpConfig->GetRelativeTolerance());
65 if (mRowForAverageOfPhiZeroed == INT_MAX)
68 this->mpLinearSystem->SetMatrixIsSymmetric(
true);
75 this->mpLinearSystem->SetMatrixIsSymmetric(
false);
88 double sqrrt_num_nodes = sqrt( this->mpMesh->GetNumNodes());
89 double normalised_basis_value = 1.0 / sqrrt_num_nodes;
99 index != dist_null_basis.
End();
102 null_basis_stripe_0[index] = 0.0;
103 null_basis_stripe_1[index] = 0.0;
104 null_basis_stripe_2[index] = normalised_basis_value;
110template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
113 if (!(this->mpBoundaryConditions->HasDirichletBoundaryConditions()))
116 if (mRowForAverageOfPhiZeroed==INT_MAX)
119 if (!mNullSpaceCreated)
122 Vec nullbasis[] = {GenerateNullBasis()};
124 this->mpLinearSystem->SetNullBasis(nullbasis, 1);
127 mNullSpaceCreated =
true;
133 this->mpLinearSystem->SetKspType(
"gmres");
134 mpConfig->SetKSPSolver(
"gmres",
true);
138 unsigned matrix_size = this->mpLinearSystem->GetSize();
139 if (!this->mMatrixIsAssembled)
143 std::vector<unsigned> row_for_average;
144 row_for_average.push_back(mRowForAverageOfPhiZeroed);
145 this->mpLinearSystem->ZeroMatrixRowsWithValueOnDiagonal(row_for_average, 0.0);
146 for (
unsigned col_index=0; col_index<matrix_size; col_index++)
148 if (((col_index-2)%3 == 0) && (col_index>1))
150 this->mpLinearSystem->SetMatrixElement(mRowForAverageOfPhiZeroed, col_index, 1);
154 this->mpLinearSystem->FinaliseLhsMatrix();
157 this->mpLinearSystem->SetRhsVectorElement(mRowForAverageOfPhiZeroed, 0);
159 this->mpLinearSystem->FinaliseRhsVector();
162 CheckCompatibilityCondition();
165template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
168 if (this->mpBoundaryConditions->HasDirichletBoundaryConditions() || mRowForAverageOfPhiZeroed!=INT_MAX )
175 DistributedVector distributed_rhs=this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(this->mpLinearSystem->rGetRhsVector());
180 index!= distributed_rhs.
End();
183 local_sum += distributed_rhs_phi_entries[index];
187 int mpi_ret = MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
188 assert(mpi_ret == MPI_SUCCESS);
190 if (fabs(global_sum)>1e-6)
194 std::cout <<
"Sum of b_{every 3 items} = " << global_sum <<
" (should be zero for compatibility)\n";
195 EXCEPTION(
"Linear system does not satisfy compatibility constraint!");
201template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
208 mBathSimulation(bathSimulation),
209 mpExtendedBidomainTissue(pTissue),
210 mpBoundaryConditions(pBcc)
212 assert(pTissue != NULL);
213 assert(pBcc != NULL);
226template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
231template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
233 std::vector<unsigned> fixedExtracellularPotentialNodes)
235 for (
unsigned i=0; i<fixedExtracellularPotentialNodes.size(); i++)
237 if (fixedExtracellularPotentialNodes[i] >= this->mpMesh->GetNumNodes() )
239 EXCEPTION(
"Fixed node number must be less than total number nodes");
243 mFixedExtracellularPotentialNodes = fixedExtracellularPotentialNodes;
246 this->mpBoundaryConditions->ResetDirichletCommunication();
248 for (
unsigned i=0; i<mFixedExtracellularPotentialNodes.size(); i++)
250 if (this->mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(mFixedExtracellularPotentialNodes[i]))
255 Node<SPACE_DIM>* p_node = this->mpMesh->GetNode(mFixedExtracellularPotentialNodes[i]);
258 GetBoundaryConditions()->AddDirichletBoundaryCondition(p_node, p_boundary_condition, 2,
false);
263template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
267 if (((row-2)%3 != 0) && row!=INT_MAX)
269 EXCEPTION(
"Row for applying the constraint 'Average of phi_e = zero' should be every 3 rows");
272 mRowForAverageOfPhiZeroed = row;
void InitialiseForSolve(Vec initialSolution)
virtual ~AbstractExtendedBidomainSolver()
void SetFixedExtracellularPotentialNodes(std::vector< unsigned > fixedExtracellularPotentialNodes)
ExtendedBidomainAssembler< ELEMENT_DIM, SPACE_DIM > * mpExtendedBidomainAssembler
virtual void FinaliseLinearSystem(Vec existingSolution)
unsigned mRowForAverageOfPhiZeroed
void PrepareForSetupLinearSystem(Vec existingSolution)
virtual void CheckCompatibilityCondition()
void SetRowForAverageOfPhiZeroed(unsigned rowMeanPhiEZero)
AbstractExtendedBidomainSolver(bool bathSimulation, AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, ExtendedBidomainTissue< SPACE_DIM > *pTissue, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, 3 > *pBcc)
virtual Vec GenerateNullBasis() const