00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "MonodomainProblem.hpp"
00030
00031 #include "Exception.hpp"
00032 #include "ReplicatableVector.hpp"
00033 #include "MonodomainDg0Assembler.hpp"
00034 #include "MonodomainMatrixBasedAssembler.hpp"
00035
00036
00037 template<unsigned DIM>
00038 AbstractCardiacPde<DIM>* MonodomainProblem<DIM>::CreateCardiacPde()
00039 {
00040 mpMonodomainPde = new MonodomainPde<DIM>(this->mpCellFactory);
00041 return mpMonodomainPde;
00042 }
00043
00044 template<unsigned DIM>
00045 AbstractDynamicAssemblerMixin<DIM, DIM, 1>* MonodomainProblem<DIM>::CreateAssembler()
00046 {
00047 assert(mpMonodomainPde);
00048
00049 if(!this->mUseMatrixBasedRhsAssembly)
00050 {
00051 MonodomainDg0Assembler<DIM,DIM>* p_assembler
00052 = new MonodomainDg0Assembler<DIM,DIM>(this->mpMesh,
00053 mpMonodomainPde,
00054 this->mpBoundaryConditionsContainer,
00055 2);
00056 return p_assembler;
00057 }
00058 else
00059 {
00060 MonodomainMatrixBasedAssembler<DIM,DIM>* p_assembler
00061 = new MonodomainMatrixBasedAssembler<DIM,DIM>(this->mpMesh,
00062 mpMonodomainPde,
00063 this->mpBoundaryConditionsContainer,
00064 2);
00065 return p_assembler;
00066 }
00067 }
00068
00069 template<unsigned DIM>
00070 MonodomainProblem<DIM>::MonodomainProblem(AbstractCardiacCellFactory<DIM>* pCellFactory)
00071 : AbstractCardiacProblem<DIM, 1>(pCellFactory),
00072 mpMonodomainPde(NULL)
00073 {
00074 }
00075
00076 template<unsigned DIM>
00077 MonodomainProblem<DIM>::~MonodomainProblem()
00078 {
00079 }
00080
00081 template<unsigned DIM>
00082 MonodomainPde<DIM> * MonodomainProblem<DIM>::GetMonodomainPde()
00083 {
00084 assert(mpMonodomainPde != NULL);
00085 return mpMonodomainPde;
00086 }
00087
00088 template<unsigned DIM>
00089 void MonodomainProblem<DIM>::WriteInfo(double time)
00090 {
00091 std::cout << "Solved to time " << time << "\n" << std::flush;
00092 ReplicatableVector voltage_replicated;
00093 voltage_replicated.ReplicatePetscVector(this->mSolution);
00094 double v_max = -DBL_MAX, v_min = DBL_MAX;
00095 for (unsigned i=0; i<this->mpMesh->GetNumNodes(); i++)
00096 {
00097 double v=voltage_replicated[i];
00098 #define COVERAGE_IGNORE
00099 if (isnan(v))
00100 {
00101 EXCEPTION("Not-a-number encountered");
00102 }
00103 #undef COVERAGE_IGNORE
00104 if ( v > v_max)
00105 {
00106 v_max = v;
00107 }
00108 if ( v < v_min)
00109 {
00110 v_min = v;
00111 }
00112 }
00113 std::cout << " V = " << "[" <<v_min << ", " << v_max << "]" << "\n" << std::flush;
00114 }
00115
00116
00118
00120
00121 template class MonodomainProblem<1>;
00122 template class MonodomainProblem<2>;
00123 template class MonodomainProblem<3>;