| Line |
Branch |
Exec |
Source |
| 1 |
|
|
// ______ ______ _ _ _____ ______ |
| 2 |
|
|
// | ____| ____| | (_)/ ____| | ____| |
| 3 |
|
|
// | |__ | |__ | | _| (___ ___| |__ |
| 4 |
|
|
// | __| | __| | | | |\___ \ / __| __| |
| 5 |
|
|
// | | | |____| |____| |____) | (__| |____ |
| 6 |
|
|
// |_| |______|______|_|_____/ \___|______| |
| 7 |
|
|
// Finite Elements for Life Sciences and Engineering |
| 8 |
|
|
// |
| 9 |
|
|
// License: LGL2.1 License |
| 10 |
|
|
// FELiScE default license: LICENSE in root folder |
| 11 |
|
|
// |
| 12 |
|
|
// Main authors: E. Schenone |
| 13 |
|
|
// |
| 14 |
|
|
|
| 15 |
|
|
// System includes |
| 16 |
|
|
|
| 17 |
|
|
// External includes |
| 18 |
|
|
|
| 19 |
|
|
// Project includes |
| 20 |
|
|
#include "Core/felisceTools.hpp" |
| 21 |
|
|
#include "Solver/linearProblemBidomainThorax.hpp" |
| 22 |
|
|
#include "InputOutput/io.hpp" |
| 23 |
|
|
#include "FiniteElement/elementMatrix.hpp" |
| 24 |
|
|
|
| 25 |
|
|
namespace felisce { |
| 26 |
|
✗ |
LinearProblemBidomainThorax::LinearProblemBidomainThorax(): |
| 27 |
|
|
LinearProblem("Problem cardiac Thorax"), |
| 28 |
|
✗ |
m_fePotThorax(nullptr), |
| 29 |
|
✗ |
m_potThoraxRes(nullptr), |
| 30 |
|
✗ |
_ECGnode(nullptr) { |
| 31 |
|
|
|
| 32 |
|
✗ |
if ( (FelisceParam::instance().writeMatrixECG) ) { |
| 33 |
|
✗ |
this->setNumberOfMatrix(2); |
| 34 |
|
|
} |
| 35 |
|
|
|
| 36 |
|
✗ |
if ( FelisceParam::instance().withCVG ) { |
| 37 |
|
|
#ifdef FELISCE_WITH_CVGRAPH |
| 38 |
|
✗ |
m_cvgDirichletVariable = "cvgraphPOTENTIAL"; |
| 39 |
|
|
#endif |
| 40 |
|
✗ |
m_seqVecs.Init("cvgraphPOTENTIAL"); |
| 41 |
|
✗ |
m_seqVecs.Init("cvgraphCURRENT"); |
| 42 |
|
|
|
| 43 |
|
|
} |
| 44 |
|
|
|
| 45 |
|
|
} |
| 46 |
|
|
|
| 47 |
|
✗ |
LinearProblemBidomainThorax::~LinearProblemBidomainThorax() { |
| 48 |
|
✗ |
m_fstransient = nullptr; |
| 49 |
|
✗ |
m_fePotThorax = nullptr; |
| 50 |
|
✗ |
delete [] m_potThoraxRes; |
| 51 |
|
✗ |
if (FelisceParam::instance().writeMatrixECG) { |
| 52 |
|
✗ |
m_thoraxRes.destroy(); |
| 53 |
|
|
} |
| 54 |
|
✗ |
if (m_boundaryConditionList.numRobinBoundaryCondition()) |
| 55 |
|
✗ |
m_robinVec.destroy(); |
| 56 |
|
✗ |
if (_ECGnode) { |
| 57 |
|
✗ |
delete [] _ECGnode; |
| 58 |
|
|
} |
| 59 |
|
|
} |
| 60 |
|
|
|
| 61 |
|
✗ |
void LinearProblemBidomainThorax::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) { |
| 62 |
|
✗ |
LinearProblem::initialize(mesh,comm, doUseSNES); |
| 63 |
|
✗ |
m_fstransient = fstransient; |
| 64 |
|
|
|
| 65 |
|
✗ |
std::vector<PhysicalVariable> listVariable(1); |
| 66 |
|
✗ |
std::vector<std::size_t> listNumComp(1); |
| 67 |
|
✗ |
listVariable[0] = potThorax; |
| 68 |
|
✗ |
listNumComp[0] = 1; |
| 69 |
|
|
//define unknown of the linear system. |
| 70 |
|
✗ |
m_listUnknown.push_back(potThorax); |
| 71 |
|
✗ |
definePhysicalVariable(listVariable,listNumComp); |
| 72 |
|
|
} |
| 73 |
|
|
|
| 74 |
|
✗ |
void LinearProblemBidomainThorax::readData(IO& io) { |
| 75 |
|
|
IGNORE_UNUSED_ARGUMENT(io); |
| 76 |
|
✗ |
m_potThoraxRes = new double[m_mesh[m_currentMesh]->numPoints()]; |
| 77 |
|
|
} |
| 78 |
|
|
|
| 79 |
|
✗ |
void LinearProblemBidomainThorax::writeSolution(int rank, std::vector<IO::Pointer>& io, double& time, int iteration) { |
| 80 |
|
✗ |
LinearProblem::writeSolution(rank,io, time, iteration); |
| 81 |
|
✗ |
if ( (FelisceParam::instance().writeMatrixECG) ) { |
| 82 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Add new variable : potThoraxRes.\n"); |
| 83 |
|
✗ |
io[m_currentMesh]->writeSolution(rank, time, iteration, 0, "potThoraxRes", m_potThoraxRes, m_supportDofUnknown[0].listNode().size(), this->dimension()); |
| 84 |
|
|
} |
| 85 |
|
|
} |
| 86 |
|
|
|
| 87 |
|
✗ |
void LinearProblemBidomainThorax::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS) { |
| 88 |
|
|
IGNORE_UNUSED_FLAG_MATRIX_RHS; |
| 89 |
|
|
IGNORE_UNUSED_ELT_TYPE; |
| 90 |
|
|
|
| 91 |
|
|
//Init pointer on Finite Element, Variable or idVariable |
| 92 |
|
✗ |
m_ipotThorax = m_listVariable.getVariableIdList(potThorax); |
| 93 |
|
✗ |
m_fePotThorax = m_listCurrentFiniteElement[m_ipotThorax]; |
| 94 |
|
|
} |
| 95 |
|
|
|
| 96 |
|
|
|
| 97 |
|
✗ |
void LinearProblemBidomainThorax::finalizeEssBCTransientDerivedProblem() { |
| 98 |
|
✗ |
std::cout << "CALL TO LinearProblemBidomainThorax::finalizeEssBCTransientDerivedProblem()"<< std::endl; |
| 99 |
|
✗ |
std::cout << "WARNING: THIS DIRECTLY READS NODEWISE DATA FROM AN ENSIGHT FILE REGARDLESS OF THE LABEL ORDERNING IN THE DATA FILE"<<std::endl; |
| 100 |
|
✗ |
bool readPotentialFromFile = false; |
| 101 |
|
|
BoundaryCondition* BC; |
| 102 |
|
✗ |
for(std::size_t iBC=0; iBC < m_boundaryConditionList.numDirichletBoundaryCondition(); iBC++) { |
| 103 |
|
✗ |
BC = m_boundaryConditionList.Dirichlet(iBC); |
| 104 |
|
✗ |
if(BC->typeValueBC() == EnsightFile) |
| 105 |
|
✗ |
readPotentialFromFile = true; |
| 106 |
|
|
} |
| 107 |
|
|
|
| 108 |
|
✗ |
if (m_boundaryConditionList.numDirichletBoundaryCondition()) { |
| 109 |
|
✗ |
PetscVector tmpRHS; |
| 110 |
|
✗ |
tmpRHS.duplicateFrom(this->sequentialSolution()); |
| 111 |
|
✗ |
tmpRHS.set(0.); |
| 112 |
|
|
|
| 113 |
|
✗ |
if (readPotentialFromFile) { |
| 114 |
|
|
felInt numBcNodeHeart; |
| 115 |
|
|
|
| 116 |
|
✗ |
if (m_thoraxMatchNode.size() == 0) { |
| 117 |
|
✗ |
m_numBcNodeThorax = readMatch(m_thoraxMatchNode,FelisceParam::instance().ECGThoraxmatchFile); |
| 118 |
|
|
} |
| 119 |
|
✗ |
if (m_heartMatchNode.size() == 0) { |
| 120 |
|
✗ |
numBcNodeHeart = readMatch(m_heartMatchNode,FelisceParam::instance().ECGmatchFile); |
| 121 |
|
✗ |
if (m_numBcNodeThorax != numBcNodeHeart) { |
| 122 |
|
✗ |
FEL_ERROR(" ERROR : Number of match nodes in Thorax different from number of match nodes in Heart.\n"); |
| 123 |
|
|
} |
| 124 |
|
|
} |
| 125 |
|
|
|
| 126 |
|
✗ |
std::vector<double> scalarValueOnDof; |
| 127 |
|
✗ |
readHeartValue(scalarValueOnDof); |
| 128 |
|
|
|
| 129 |
|
|
felInt iPos; |
| 130 |
|
|
double value; |
| 131 |
|
✗ |
felInt sizeVent = 45580; |
| 132 |
|
|
|
| 133 |
|
✗ |
for (felInt i = 0; i < m_numBcNodeThorax; i++) { |
| 134 |
|
✗ |
iPos = m_thoraxMatchNode[i]-1; |
| 135 |
|
✗ |
AOApplicationToPetsc(m_ao,1,&iPos); |
| 136 |
|
✗ |
double penalty = FelisceParam::instance().alphaRobin[0]; // ventricles |
| 137 |
|
✗ |
if (m_heartMatchNode[i]-1 > (sizeVent-1) ) { |
| 138 |
|
✗ |
penalty = FelisceParam::instance().alphaRobin[1]; // atria |
| 139 |
|
|
} |
| 140 |
|
✗ |
value = scalarValueOnDof[m_heartMatchNode[i]-1]*penalty; |
| 141 |
|
✗ |
tmpRHS.setValue(iPos,value,INSERT_VALUES); |
| 142 |
|
|
} |
| 143 |
|
✗ |
tmpRHS.assembly(); |
| 144 |
|
|
} |
| 145 |
|
|
|
| 146 |
|
|
// Use the heart extra-cell potential solution to enforce the Dirichlet |
| 147 |
|
|
double solval; |
| 148 |
|
✗ |
std::map<int, double> idBCAndValue; |
| 149 |
|
✗ |
std::set<felInt> idDofBC; |
| 150 |
|
|
|
| 151 |
|
|
int rankProc; |
| 152 |
|
✗ |
MPI_Comm_rank(PETSC_COMM_WORLD, &rankProc); |
| 153 |
|
|
|
| 154 |
|
|
//BoundaryCondition* BC; |
| 155 |
|
✗ |
for(std::size_t iBC=0; iBC < m_boundaryConditionList.numDirichletBoundaryCondition(); iBC++) { |
| 156 |
|
✗ |
BC = m_boundaryConditionList.Dirichlet(iBC); |
| 157 |
|
|
|
| 158 |
|
✗ |
if(BC->typeValueBC() == EnsightFile) { |
| 159 |
|
|
|
| 160 |
|
✗ |
felInt idEltInSupportLocal = 0; |
| 161 |
|
✗ |
felInt idEltInSupportGlobal = 0; |
| 162 |
|
|
felInt idDof; |
| 163 |
|
✗ |
std::set<felInt> idDof2; |
| 164 |
|
|
|
| 165 |
|
✗ |
std::vector< std::pair<felInt,felInt> > idEltAndIdSupport; |
| 166 |
|
✗ |
idEltAndIdSupport = BC->idEltAndIdSupport(); |
| 167 |
|
|
|
| 168 |
|
✗ |
std::vector<double>& valueBCInSupportDof = BC->ValueBCInSupportDof(); |
| 169 |
|
✗ |
valueBCInSupportDof.resize(idEltAndIdSupport.size(), 0.); |
| 170 |
|
|
|
| 171 |
|
✗ |
for (std::size_t iSupportDofBC = 0; iSupportDofBC < idEltAndIdSupport.size(); iSupportDofBC++) { |
| 172 |
|
✗ |
idEltInSupportLocal = idEltAndIdSupport[iSupportDofBC].first; |
| 173 |
|
✗ |
ISLocalToGlobalMappingApply(m_mappingElem[m_currentMesh], 1, &idEltInSupportLocal, &idEltInSupportGlobal); |
| 174 |
|
✗ |
dof().loc2glob(idEltInSupportGlobal, idEltAndIdSupport[iSupportDofBC].second, 0, 0, idDof); |
| 175 |
|
✗ |
AOApplicationToPetsc(m_ao, 1, &idDof); |
| 176 |
|
✗ |
bool is_new_element = idDof2.insert(idDof).second; |
| 177 |
|
|
|
| 178 |
|
✗ |
if (is_new_element) { |
| 179 |
|
✗ |
tmpRHS.getValues(1, &idDof, &solval); |
| 180 |
|
✗ |
valueBCInSupportDof[iSupportDofBC] = solval; |
| 181 |
|
|
} |
| 182 |
|
|
} |
| 183 |
|
|
} |
| 184 |
|
|
} |
| 185 |
|
✗ |
tmpRHS.destroy(); |
| 186 |
|
|
} |
| 187 |
|
|
} |
| 188 |
|
|
|
| 189 |
|
✗ |
void LinearProblemBidomainThorax::allocateVectorBoundaryConditionDerivedLinPb() { |
| 190 |
|
✗ |
bool readPotentialFromFile = false; |
| 191 |
|
|
BoundaryCondition* BC; |
| 192 |
|
✗ |
for(std::size_t iBC=0; iBC < m_boundaryConditionList.numRobinBoundaryCondition(); iBC++) { |
| 193 |
|
✗ |
BC = m_boundaryConditionList.Robin(iBC); |
| 194 |
|
✗ |
if(BC->typeValueBC() == EnsightFile) |
| 195 |
|
✗ |
readPotentialFromFile = true; |
| 196 |
|
|
} |
| 197 |
|
|
|
| 198 |
|
✗ |
if (readPotentialFromFile) { |
| 199 |
|
|
felInt numBcNodeHeart; |
| 200 |
|
|
|
| 201 |
|
✗ |
if (m_thoraxMatchNode.size() == 0) { |
| 202 |
|
✗ |
m_numBcNodeThorax = readMatch(m_thoraxMatchNode,FelisceParam::instance().ECGThoraxmatchFile); |
| 203 |
|
|
} |
| 204 |
|
✗ |
if (m_heartMatchNode.size() == 0) { |
| 205 |
|
✗ |
numBcNodeHeart = readMatch(m_heartMatchNode,FelisceParam::instance().ECGmatchFile); |
| 206 |
|
✗ |
if (m_numBcNodeThorax != numBcNodeHeart) { |
| 207 |
|
✗ |
FEL_ERROR(" ERROR : Number of match nodes in Thorax different from number of match nodes in Heart.\n"); |
| 208 |
|
|
} |
| 209 |
|
|
} |
| 210 |
|
|
|
| 211 |
|
✗ |
std::vector<double> scalarValueOnDof; |
| 212 |
|
✗ |
readHeartValue(scalarValueOnDof); |
| 213 |
|
|
|
| 214 |
|
✗ |
if (m_fstransient->iteration == 1) { |
| 215 |
|
✗ |
m_robinVec.duplicateFrom(this->sequentialSolution()); |
| 216 |
|
|
} |
| 217 |
|
|
|
| 218 |
|
✗ |
m_robinVec.set(0.); |
| 219 |
|
|
|
| 220 |
|
|
felInt iPos; |
| 221 |
|
|
double value; |
| 222 |
|
|
|
| 223 |
|
✗ |
for (felInt i = 0; i < m_numBcNodeThorax; i++) { |
| 224 |
|
✗ |
iPos = m_thoraxMatchNode[i]-1; |
| 225 |
|
✗ |
AOApplicationToPetsc(m_ao,1,&iPos); |
| 226 |
|
✗ |
value = scalarValueOnDof[m_heartMatchNode[i]-1]; |
| 227 |
|
✗ |
m_robinVec.setValue(iPos,value,INSERT_VALUES); |
| 228 |
|
|
} |
| 229 |
|
✗ |
m_robinVec.assembly(); |
| 230 |
|
|
} |
| 231 |
|
|
} |
| 232 |
|
|
|
| 233 |
|
✗ |
void LinearProblemBidomainThorax::readElectrodesNode(bool flagTranspos) { |
| 234 |
|
|
/* Leads position rules: |
| 235 |
|
|
|
| 236 |
|
|
V1 est placée sur le 4ème espace intercostal droit, au bord droit du sternum. |
| 237 |
|
|
V2 est placée sur le 4ème espace intercostal gauche, au bord gauche du sternum. |
| 238 |
|
|
V4 est placée sur le 5ème espace intercostal gauche, sur la ligne médioclaviculaire. |
| 239 |
|
|
V3 est placée entre V2 et V4. |
| 240 |
|
|
V5 est placée sur le 5ème espace intercostal gauche, sur la ligne axillaire antérieure. |
| 241 |
|
|
V6 est placée sur le 5ème espace intercostal gauche, sur la ligne axillaire moyenne. |
| 242 |
|
|
*/ |
| 243 |
|
|
|
| 244 |
|
|
// Standard transpos transfer matrix |
| 245 |
|
✗ |
if (flagTranspos) { |
| 246 |
|
✗ |
std::string electrodeFileName; |
| 247 |
|
✗ |
electrodeFileName = FelisceParam::instance().inputDirectory + "/" + FelisceParam::instance().ECGelectrodeFile;// load the electrode file (thorax skin nodes) |
| 248 |
|
|
|
| 249 |
|
✗ |
std::ifstream electrodeFile(electrodeFileName.c_str()); |
| 250 |
|
✗ |
if(!electrodeFile) { |
| 251 |
|
✗ |
FEL_ERROR("Fatal error: cannot open match file " + electrodeFileName ); |
| 252 |
|
|
} |
| 253 |
|
✗ |
electrodeFile >> m_numSourceNode; // For standard 12 leads ECG = 9 |
| 254 |
|
|
|
| 255 |
|
✗ |
_ECGnode = new felInt[m_numSourceNode]; |
| 256 |
|
|
|
| 257 |
|
✗ |
for(felInt j=0; j<m_numSourceNode; j++) { |
| 258 |
|
✗ |
electrodeFile >> _ECGnode[j]; |
| 259 |
|
|
} |
| 260 |
|
✗ |
electrodeFile.close(); |
| 261 |
|
|
} |
| 262 |
|
|
// To calculate the transfer matrix on the whole heart surface, reading match files |
| 263 |
|
|
else { |
| 264 |
|
✗ |
std::string matchFileName; |
| 265 |
|
✗ |
matchFileName = FelisceParam::instance().inputDirectory + FelisceParam::instance().ECGmatchFile;// load the match file (heart nodes) |
| 266 |
|
|
|
| 267 |
|
✗ |
std::ifstream matchFile(matchFileName.c_str()); |
| 268 |
|
✗ |
if(!matchFile) { |
| 269 |
|
✗ |
FEL_ERROR("Fatal error: cannot open match file " + matchFileName ); |
| 270 |
|
|
} |
| 271 |
|
✗ |
matchFile >> m_numSourceNode; |
| 272 |
|
|
|
| 273 |
|
✗ |
_ECGnode = new felInt[m_numSourceNode]; |
| 274 |
|
|
|
| 275 |
|
✗ |
for(felInt j=0; j<m_numSourceNode; j++) { |
| 276 |
|
✗ |
matchFile >> _ECGnode[j]; |
| 277 |
|
|
} |
| 278 |
|
✗ |
matchFile.close(); |
| 279 |
|
|
} |
| 280 |
|
|
|
| 281 |
|
|
} |
| 282 |
|
|
|
| 283 |
|
✗ |
void LinearProblemBidomainThorax::calculateRHS() { |
| 284 |
|
✗ |
if (m_fstransient->iteration == 1) { |
| 285 |
|
✗ |
readElectrodesNode(); |
| 286 |
|
|
} |
| 287 |
|
|
|
| 288 |
|
✗ |
double itApply = m_fstransient->time; |
| 289 |
|
|
felInt iPos; |
| 290 |
|
|
|
| 291 |
|
✗ |
vector().set(0.); |
| 292 |
|
|
// Impose a delta_x(x_i) for i = 1, ..., numSourceNode (at time = i) |
| 293 |
|
✗ |
if ( ( static_cast<int>(itApply) <= m_numSourceNode ) && ( Tools::equal(fmod(itApply,1),0.) ) ) { |
| 294 |
|
✗ |
int i = _ECGnode[static_cast<int>(itApply)-1]-1; |
| 295 |
|
✗ |
iPos = i; |
| 296 |
|
✗ |
AOPetscToApplication(m_ao,1,&iPos); |
| 297 |
|
✗ |
vector().setValue(iPos,1.0,INSERT_VALUES); |
| 298 |
|
|
} |
| 299 |
|
|
|
| 300 |
|
✗ |
vector().assembly(); |
| 301 |
|
|
|
| 302 |
|
|
} |
| 303 |
|
|
|
| 304 |
|
✗ |
void LinearProblemBidomainThorax::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) { |
| 305 |
|
|
IGNORE_UNUSED_IEL; |
| 306 |
|
|
IGNORE_UNUSED_FLAG_MATRIX_RHS; |
| 307 |
|
|
IGNORE_UNUSED_ELEM_ID_POINT; |
| 308 |
|
|
|
| 309 |
|
✗ |
m_fePotThorax->updateMeasQuadPt(0, elemPoint); |
| 310 |
|
✗ |
m_fePotThorax->updateFirstDeriv(0, elemPoint); |
| 311 |
|
|
|
| 312 |
|
|
// Assembling matrix(0) |
| 313 |
|
|
// Evaluate conductivity. |
| 314 |
|
|
double sigmaTval; |
| 315 |
|
✗ |
ThoraxConductivity* fctSigmaT = new ThoraxConductivity(); |
| 316 |
|
✗ |
fctSigmaT->initialize(m_fstransient); |
| 317 |
|
✗ |
evalFunctionOnElem(*fctSigmaT,sigmaTval,m_currentLabel); |
| 318 |
|
|
// \sigma_T(x) \grad u_T |
| 319 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(sigmaTval,*m_fePotThorax,0,0,1); |
| 320 |
|
✗ |
if ( (FelisceParam::instance().writeMatrixECG) ) { |
| 321 |
|
|
// matrix(1) = matrix(0) without BC (used to compute residual) |
| 322 |
|
✗ |
m_elementMat[1]->grad_phi_i_grad_phi_j(sigmaTval,*m_fePotThorax,0,0,1); |
| 323 |
|
|
} |
| 324 |
|
✗ |
delete fctSigmaT; |
| 325 |
|
|
|
| 326 |
|
|
} |
| 327 |
|
|
|
| 328 |
|
✗ |
void LinearProblemBidomainThorax::userElementComputeNaturalBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>&, felInt& iel, int) { |
| 329 |
|
|
//std::cout << "LinearProblemBidomainThorax::userElementComputeNaturalBoundaryCondition" << std::endl; |
| 330 |
|
✗ |
if (m_boundaryConditionList.numRobinBoundaryCondition()) { |
| 331 |
|
|
BoundaryCondition* BC; |
| 332 |
|
✗ |
CurvilinearFiniteElement* fe = m_listCurvilinearFiniteElement[0]; |
| 333 |
|
✗ |
fe -> updateMeas(0, elemPoint); |
| 334 |
|
|
|
| 335 |
|
✗ |
for (std::size_t iRobin=0; iRobin<m_boundaryConditionList.numRobinBoundaryCondition(); iRobin++ ) { |
| 336 |
|
✗ |
if ( m_elemFieldRobin[iRobin].getType() != CONSTANT_FIELD ) { |
| 337 |
|
✗ |
BC = m_boundaryConditionList.Robin(iRobin); |
| 338 |
|
✗ |
for(auto it_labelNumber = BC->listLabel().begin(); it_labelNumber != BC->listLabel().end(); it_labelNumber++) { |
| 339 |
|
✗ |
if(*it_labelNumber == m_currentLabel) { |
| 340 |
|
✗ |
m_elemFieldRobin[iRobin].setValue(m_robinVec, *fe, iel, 0, m_ao, dof()); |
| 341 |
|
✗ |
double coef = FelisceParam::instance().alphaRobin[iRobin]; |
| 342 |
|
✗ |
FelisceParam::instance().betaRobin[iRobin] = 1./coef; |
| 343 |
|
|
//m_elementVectorBD[0]->sourceForComp(coef,*fe,m_elemFieldRobin[iRobin],0,0); |
| 344 |
|
|
} |
| 345 |
|
|
} |
| 346 |
|
|
} |
| 347 |
|
|
} |
| 348 |
|
|
} |
| 349 |
|
|
|
| 350 |
|
|
} |
| 351 |
|
|
|
| 352 |
|
✗ |
felInt LinearProblemBidomainThorax::readMatch(std::vector<felInt>& matchNode, std::string matchName) { |
| 353 |
|
|
|
| 354 |
|
|
felInt numSourceNode; |
| 355 |
|
|
// Input files. |
| 356 |
|
✗ |
std::string matchFileName = FelisceParam::instance().inputDirectory + "/" + matchName;// load the match file (heart nodes) |
| 357 |
|
|
|
| 358 |
|
✗ |
std::ifstream matchFile(matchFileName.c_str()); |
| 359 |
|
✗ |
if(!matchFile) { |
| 360 |
|
✗ |
FEL_ERROR("Fatal error: cannot open match file " + matchFileName ); |
| 361 |
|
|
} |
| 362 |
|
✗ |
if (FelisceParam::verbose() > 0) { |
| 363 |
|
✗ |
std::cout << "Reading " << matchFileName << std::endl; |
| 364 |
|
|
} |
| 365 |
|
✗ |
matchFile >> numSourceNode; |
| 366 |
|
✗ |
matchNode.resize(numSourceNode); |
| 367 |
|
|
felInt value; |
| 368 |
|
✗ |
for(int j=0; j<numSourceNode; j++) { |
| 369 |
|
✗ |
matchFile >> value; |
| 370 |
|
✗ |
matchNode[j]=value; |
| 371 |
|
|
} |
| 372 |
|
✗ |
matchFile.close(); |
| 373 |
|
|
|
| 374 |
|
✗ |
return numSourceNode; |
| 375 |
|
|
|
| 376 |
|
|
} |
| 377 |
|
|
|
| 378 |
|
✗ |
void LinearProblemBidomainThorax::readHeartValue(std::vector<double>& scalarValue) { |
| 379 |
|
|
|
| 380 |
|
✗ |
std::ostringstream index; |
| 381 |
|
✗ |
index << m_fstransient->iteration; |
| 382 |
|
✗ |
std::string strindex = ""; |
| 383 |
|
✗ |
strindex = index.str(); |
| 384 |
|
✗ |
std::string varName = "potExtraCell"; |
| 385 |
|
|
// std::string varName = "potThorax"; |
| 386 |
|
✗ |
std::string sclFile; |
| 387 |
|
|
|
| 388 |
|
✗ |
if ( m_fstransient->iteration < 10 ) { |
| 389 |
|
✗ |
sclFile = varName + ".0000" + strindex; |
| 390 |
|
✗ |
} else if ( m_fstransient->iteration < 100 ) { |
| 391 |
|
✗ |
sclFile = varName + ".000" + strindex; |
| 392 |
|
✗ |
} else if ( m_fstransient->iteration < 1000 ) { |
| 393 |
|
✗ |
sclFile = varName + ".00" + strindex; |
| 394 |
|
✗ |
} else if ( m_fstransient->iteration < 10000 ) { |
| 395 |
|
✗ |
sclFile = varName + ".0" + strindex; |
| 396 |
|
|
} else { |
| 397 |
|
✗ |
sclFile = varName + "." + strindex; |
| 398 |
|
|
} |
| 399 |
|
✗ |
std::string bcFileName = FelisceParam::instance().bcCondDir + "/" + sclFile + ".scl"; |
| 400 |
|
✗ |
std::ifstream bcFile(bcFileName.c_str(), std::ios_base::in); |
| 401 |
|
|
|
| 402 |
|
✗ |
if ( !bcFile ) { |
| 403 |
|
✗ |
FEL_ERROR(" ERROR: Can not open file "+bcFileName+"."); |
| 404 |
|
|
} else { |
| 405 |
|
✗ |
if (FelisceParam::verbose() > 0) { |
| 406 |
|
✗ |
std::cout << "Reading " << bcFileName << std::endl; |
| 407 |
|
|
} |
| 408 |
|
|
} |
| 409 |
|
|
|
| 410 |
|
|
char newline[1024]; |
| 411 |
|
✗ |
bcFile.getline(newline, 1024, '\n'); |
| 412 |
|
|
double value; |
| 413 |
|
✗ |
felInt i = 0; |
| 414 |
|
✗ |
while (!bcFile.eof()) { |
| 415 |
|
✗ |
bcFile >> value; |
| 416 |
|
✗ |
scalarValue.push_back(value); |
| 417 |
|
✗ |
i++; |
| 418 |
|
|
} |
| 419 |
|
✗ |
bcFile.close(); |
| 420 |
|
|
} |
| 421 |
|
|
|
| 422 |
|
✗ |
void LinearProblemBidomainThorax::calculateRes() { |
| 423 |
|
✗ |
if (m_fstransient->iteration == 1) { |
| 424 |
|
✗ |
m_thoraxRes.duplicateFrom(vector()); |
| 425 |
|
|
} |
| 426 |
|
✗ |
m_thoraxRes.set(0.); |
| 427 |
|
|
|
| 428 |
|
|
// m_thoraxRes = _RHS_aux - matrix(1)*m_sol |
| 429 |
|
|
|
| 430 |
|
|
// 1. m_thoraxRes = matrix(1)*this->solution() |
| 431 |
|
✗ |
mult(matrix(1),this->solution(),m_thoraxRes); |
| 432 |
|
|
|
| 433 |
|
|
// 2. m_thoraxRes = _RHS_aux - m_thoraxRes |
| 434 |
|
✗ |
double a = -1.0; |
| 435 |
|
✗ |
m_thoraxRes.scale(a); // _RHS_aux = zeros (_RHS without BC) |
| 436 |
|
|
|
| 437 |
|
✗ |
felInt* iPos = new felInt[m_numDof]; |
| 438 |
|
✗ |
for (felInt i = 0; i < m_numDof; i++) { |
| 439 |
|
✗ |
iPos[i] = i; |
| 440 |
|
|
} |
| 441 |
|
✗ |
AOPetscToApplication(m_ao,m_numDof,iPos); |
| 442 |
|
✗ |
m_thoraxRes.getValues(m_numDof,iPos,m_potThoraxRes); |
| 443 |
|
✗ |
delete [] iPos; |
| 444 |
|
|
|
| 445 |
|
|
} |
| 446 |
|
|
|
| 447 |
|
✗ |
void LinearProblemBidomainThorax::computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, felInt& iel, FlagMatrixRHS /*flagMatrixRHS*/) { |
| 448 |
|
|
|
| 449 |
|
|
(void)iel; |
| 450 |
|
✗ |
m_curvFe->updateMeas(0,elemPoint); |
| 451 |
|
|
} |
| 452 |
|
|
|
| 453 |
|
|
#ifdef FELISCE_WITH_CVGRAPH |
| 454 |
|
✗ |
void LinearProblemBidomainThorax::initMassBoundaryForCVG() { |
| 455 |
|
✗ |
this->assembleMassBoundaryOnly(&LinearProblemBidomainThorax::massMatrixComputer, |
| 456 |
|
✗ |
this->slave()->interfaceLabels(/*iConn*/0), |
| 457 |
|
|
&LinearProblemBidomainThorax::initPerETMass, |
| 458 |
|
|
&LinearProblemBidomainThorax::updateFE, |
| 459 |
|
✗ |
this->dofBD(/*iConn*/0), this->massBDVec()[/*iConn*/0]); |
| 460 |
|
|
} |
| 461 |
|
|
|
| 462 |
|
✗ |
void LinearProblemBidomainThorax::massMatrixComputer(felInt ielSupportDof) { |
| 463 |
|
✗ |
this->m_elementMatBD[0]->zero(); |
| 464 |
|
✗ |
this->m_elementMatBD[0]->phi_i_phi_j(/*coef*/1.,*m_curvFe,0,0,1); |
| 465 |
|
✗ |
this->setValueMatrixBD(ielSupportDof); |
| 466 |
|
|
} |
| 467 |
|
|
|
| 468 |
|
|
void |
| 469 |
|
✗ |
LinearProblemBidomainThorax::initPerETMass() { |
| 470 |
|
✗ |
felInt numDofTotal = 0; |
| 471 |
|
|
//pay attention this numDofTotal is bigger than expected since it counts for all the VARIABLES |
| 472 |
|
✗ |
for (std::size_t iFe = 0; iFe < this->m_listCurvilinearFiniteElement.size(); iFe++) {//this loop is on variables while it should be on unknown |
| 473 |
|
✗ |
numDofTotal += this->m_listCurvilinearFiniteElement[iFe]->numDof()*this->m_listVariable[iFe].numComponent(); |
| 474 |
|
|
} |
| 475 |
|
✗ |
m_globPosColumn.resize(numDofTotal); m_globPosRow.resize(numDofTotal); m_matrixValues.resize(numDofTotal*numDofTotal); |
| 476 |
|
|
//std::cout << "LinearProblemBidomainThorax::initPerETMass() m_ipotThorax = " << m_ipotThorax << std::endl; |
| 477 |
|
✗ |
m_curvFe = this->m_listCurvilinearFiniteElement[ m_ipotThorax ]; |
| 478 |
|
|
} |
| 479 |
|
|
|
| 480 |
|
|
void |
| 481 |
|
✗ |
LinearProblemBidomainThorax::updateFE(const std::vector<Point*>& elemPoint, const std::vector<int>&) { |
| 482 |
|
✗ |
m_curvFe->updateMeasNormal(0, elemPoint); |
| 483 |
|
|
} |
| 484 |
|
|
|
| 485 |
|
✗ |
void LinearProblemBidomainThorax::readData() { |
| 486 |
|
|
//std::cout << "LinearProblemBidomainThorax::readData()" <<std::endl; |
| 487 |
|
✗ |
for ( std::size_t iConn(0); iConn<this->slave()->numConnections(); ++iConn ) { |
| 488 |
|
✗ |
std::vector<PetscVector> vecs; |
| 489 |
|
✗ |
for ( std::size_t cVar(0); cVar<this->slave()->numVarToRead(iConn); ++cVar) { |
| 490 |
|
✗ |
std::string varToRead = this->slave()->readVariable(iConn,cVar); |
| 491 |
|
|
//std::cout << "iConn = " << iConn << " cVar = " << cVar << " varToRead = " << varToRead << std::endl; |
| 492 |
|
✗ |
if ( varToRead == "POTENTIAL" ) { |
| 493 |
|
✗ |
vecs.push_back(m_seqVecs.Get("cvgraphPOTENTIAL")); |
| 494 |
|
✗ |
} else if ( varToRead == "CURRENT" ){ |
| 495 |
|
✗ |
vecs.push_back(m_seqVecs.Get("cvgraphCURRENT")); |
| 496 |
|
|
} else { |
| 497 |
|
✗ |
std::cout<<"Connection: "<<iConn<<"--variable to read: "<< varToRead << std::endl; |
| 498 |
|
✗ |
FEL_ERROR("Error in variable to read"); |
| 499 |
|
|
} |
| 500 |
|
|
} |
| 501 |
|
✗ |
this->slave()->receiveData(vecs,iConn); |
| 502 |
|
|
} |
| 503 |
|
|
} |
| 504 |
|
|
|
| 505 |
|
✗ |
void LinearProblemBidomainThorax::sendData() { |
| 506 |
|
|
//std::cout << "LinearProblemBidomainThorax::sendData()" <<std::endl; |
| 507 |
|
✗ |
this->gatherSolution(); |
| 508 |
|
✗ |
for ( std::size_t iConn(0); iConn<this->slave()->numConnections(); ++iConn ) { |
| 509 |
|
✗ |
std::vector<PetscVector> vecs; |
| 510 |
|
✗ |
for ( std::size_t cVar(0); cVar<this->slave()->numVarToSend(iConn); ++cVar) { |
| 511 |
|
✗ |
std::string varToSend = this->slave()->sendVariable(iConn, cVar); |
| 512 |
|
✗ |
if ( varToSend == "POTENTIAL" ) { |
| 513 |
|
✗ |
vecs.push_back(this->sequentialSolution()); |
| 514 |
|
✗ |
} else if ( varToSend == "CURRENT" ) { |
| 515 |
|
✗ |
this->computeResidual(); |
| 516 |
|
✗ |
vecs.push_back(this->seqResidual()); |
| 517 |
|
|
} |
| 518 |
|
|
} |
| 519 |
|
✗ |
this->slave()->sendData(vecs,iConn); |
| 520 |
|
|
} |
| 521 |
|
|
} |
| 522 |
|
|
|
| 523 |
|
✗ |
void LinearProblemBidomainThorax::initMatrix() { |
| 524 |
|
✗ |
this->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs); |
| 525 |
|
✗ |
this->assembleMatrixRHS(MpiInfo::rankProc(),FlagMatrixRHS::matrix_and_rhs); |
| 526 |
|
✗ |
this->createAndCopyMatrixRHSWithoutBC(); |
| 527 |
|
|
} |
| 528 |
|
|
|
| 529 |
|
✗ |
void LinearProblemBidomainThorax::addCurrentToRHS() { |
| 530 |
|
|
//! for each dof on the boundary on this proc. |
| 531 |
|
✗ |
felInt numLocalDofInterface=this->dofBD(/*iBD*/0).numLocalDofInterface(); |
| 532 |
|
✗ |
std::vector<double> tmp(numLocalDofInterface); |
| 533 |
|
✗ |
m_seqVecs.Get("cvgraphCURRENT").getValues( numLocalDofInterface, this->dofBD(/*iBD*/0).loc2PetscVolPtr(), tmp.data() ); |
| 534 |
|
✗ |
this->vector().setValues(numLocalDofInterface,this->dofBD(/*iBD*/0).loc2PetscVolPtr(),tmp.data(),ADD_VALUES); |
| 535 |
|
✗ |
this->vector().assembly(); |
| 536 |
|
|
} |
| 537 |
|
|
|
| 538 |
|
✗ |
void LinearProblemBidomainThorax::cvgAddRobinToRHS(){ |
| 539 |
|
|
//std::cout << "LinearProblemBidomain::cvgAddRobinToRHS()" << std::endl; |
| 540 |
|
|
//std::cout << " alpha_robin = " << FelisceParam::instance().alphaRobin[0/*iRobin*/] << std::endl; |
| 541 |
|
|
} |
| 542 |
|
|
|
| 543 |
|
|
#endif |
| 544 |
|
|
|
| 545 |
|
|
} |
| 546 |
|
|
|