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: A. Collin |
13 |
|
|
// |
14 |
|
|
|
15 |
|
|
// System includes |
16 |
|
|
#include <iostream> |
17 |
|
|
#include <fstream> |
18 |
|
|
|
19 |
|
|
// External includes |
20 |
|
|
|
21 |
|
|
// Project includes |
22 |
|
|
#include "Solver/linearProblemBidomain.hpp" |
23 |
|
|
#include "InputOutput/io.hpp" |
24 |
|
|
#include "Core/felisceTools.hpp" |
25 |
|
|
#include "Core/filesystemUtil.hpp" |
26 |
|
|
#include "Tools/fe_utilities.hpp" |
27 |
|
|
|
28 |
|
|
namespace felisce |
29 |
|
|
{ |
30 |
|
✗ |
LinearProblemBidomain::LinearProblemBidomain(): |
31 |
|
|
LinearProblem("Problem cardiac Bidomain",2), |
32 |
|
✗ |
m_rom(nullptr), |
33 |
|
✗ |
m_fePotTransMemb(nullptr), |
34 |
|
✗ |
m_fePotExtraCell(nullptr), |
35 |
|
✗ |
m_fePotTransMembCurv(nullptr), |
36 |
|
✗ |
m_fePotExtraCellCurv(nullptr), |
37 |
|
✗ |
m_sortedSol(nullptr), |
38 |
|
✗ |
m_bdf(nullptr), |
39 |
|
✗ |
m_initializeCVGraphInterface(false), |
40 |
|
✗ |
m_electrodeNode(nullptr), |
41 |
|
✗ |
m_electrodeMeasure(nullptr) { |
42 |
|
✗ |
if (FelisceParam::instance().stateFilter) |
43 |
|
✗ |
this->setNumberOfMatrix(3); |
44 |
|
|
|
45 |
|
|
|
46 |
|
✗ |
if ( FelisceParam::instance().withCVG ) { |
47 |
|
|
#ifdef FELISCE_WITH_CVGRAPH |
48 |
|
✗ |
m_cvgDirichletVariable = "cvgraphPOTENTIAL"; |
49 |
|
|
#endif |
50 |
|
✗ |
m_seqVecs.Init("cvgraphPOTENTIAL"); |
51 |
|
✗ |
m_seqVecs.Init("cvgraphCURRENT"); |
52 |
|
|
} |
53 |
|
|
|
54 |
|
|
} |
55 |
|
|
|
56 |
|
✗ |
LinearProblemBidomain::~LinearProblemBidomain() { |
57 |
|
✗ |
m_rom = nullptr; |
58 |
|
✗ |
m_fstransient = nullptr; |
59 |
|
✗ |
m_fePotTransMemb = nullptr; |
60 |
|
✗ |
m_fePotExtraCell = nullptr; |
61 |
|
✗ |
m_bdf = nullptr; |
62 |
|
✗ |
delete [] m_sortedSol; |
63 |
|
|
// Elisa, Dec. 2014 - Luenberger filter |
64 |
|
✗ |
if (FelisceParam::instance().electrodeControl) { |
65 |
|
✗ |
delete [] m_electrodeNode; |
66 |
|
✗ |
for (felInt i=0; i<m_numElectrodeNode+1; i++) { |
67 |
|
✗ |
delete [] m_electrodeMeasure[i]; |
68 |
|
|
} |
69 |
|
✗ |
delete [] m_electrodeMeasure; |
70 |
|
✗ |
m_electrodeControl.destroy(); |
71 |
|
|
} |
72 |
|
|
} |
73 |
|
|
|
74 |
|
✗ |
void LinearProblemBidomain::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) { |
75 |
|
✗ |
LinearProblem::initialize(mesh,comm, doUseSNES); |
76 |
|
✗ |
m_fstransient = fstransient; |
77 |
|
|
|
78 |
|
✗ |
std::vector<PhysicalVariable> listVariable(2); |
79 |
|
✗ |
std::vector<std::size_t> listNumComp(2); |
80 |
|
|
|
81 |
|
✗ |
listVariable[0] = potTransMemb; |
82 |
|
✗ |
listNumComp[0] = 1; |
83 |
|
✗ |
listVariable[1] = potExtraCell; |
84 |
|
✗ |
listNumComp[1] = 1; |
85 |
|
|
//define unknown of the linear system. |
86 |
|
✗ |
m_listUnknown.push_back(potTransMemb); |
87 |
|
✗ |
m_listUnknown.push_back(potExtraCell); |
88 |
|
✗ |
definePhysicalVariable(listVariable,listNumComp); |
89 |
|
|
} |
90 |
|
|
|
91 |
|
✗ |
void LinearProblemBidomain::readData(IO& io) { |
92 |
|
|
// Read fibers file (*.00000.vct) |
93 |
|
✗ |
m_vectorFiber.resize(m_mesh[m_currentMesh]->numPoints()*3); |
94 |
|
✗ |
io.readVariable(0, 0.,&m_vectorFiber[0], m_vectorFiber.size()); |
95 |
|
|
// Read distance from endocardium file for Zygote geometry (*.00000.scl) |
96 |
|
✗ |
if ( (FelisceParam::instance().typeOfAppliedCurrent == "zygote") || (FelisceParam::instance().typeOfAppliedCurrent == "zygoteImproved") || (FelisceParam::instance().typeOfAppliedCurrent == "heart") || (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart") ) { |
97 |
|
✗ |
m_vectorDistance.resize(m_numDof,0.); |
98 |
|
✗ |
io.readVariable(1, 0.,&m_vectorDistance[0],m_mesh[m_currentMesh]->numPoints()); |
99 |
|
✗ |
if ( (FelisceParam::instance().typeOfAppliedCurrent == "heart") || (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart") ) { |
100 |
|
✗ |
if (FelisceParam::instance().hasCoupledAtriaVent) { |
101 |
|
✗ |
m_vectorAngle.resize(m_numDof,0.); |
102 |
|
✗ |
io.readVariable(2, 0.,&m_vectorAngle[0],m_mesh[m_currentMesh]->numPoints()); |
103 |
|
✗ |
m_vectorRef.resize(m_numDof,0.); |
104 |
|
✗ |
io.readVariable(3, 0.,&m_vectorRef[0],m_mesh[m_currentMesh]->numPoints()); |
105 |
|
|
} |
106 |
|
|
} |
107 |
|
|
} |
108 |
|
|
} |
109 |
|
|
|
110 |
|
✗ |
void LinearProblemBidomain::updateFibers(IO& io,double iteration) { //when the mesh moves |
111 |
|
✗ |
m_vectorFiber.resize(m_mesh[m_currentMesh]->numPoints()*3); |
112 |
|
✗ |
io.readVariable(0, iteration,&m_vectorFiber[0], m_vectorFiber.size()); |
113 |
|
|
} |
114 |
|
|
|
115 |
|
✗ |
void LinearProblemBidomain::readDataDisplacement(std::vector<IO::Pointer>& io,double iteration) { |
116 |
|
|
// Read displacement file (*.00000.vct) |
117 |
|
✗ |
if ((FelisceParam::instance().typeOfAppliedCurrent == "zygote")&&(FelisceParam::instance().typeOfAppliedCurrent == "planarWave")){ |
118 |
|
✗ |
const int iMesh = m_listVariable[2].idMesh(); |
119 |
|
✗ |
m_vectorDisp.resize(m_mesh[iMesh]->numPoints()*3,0.); |
120 |
|
✗ |
io[iMesh]->readVariable(2, iteration,&m_vectorDisp[0], m_vectorDisp.size()); |
121 |
|
|
} |
122 |
|
|
else{ |
123 |
|
✗ |
const int iMesh = m_listVariable[1].idMesh(); |
124 |
|
✗ |
m_vectorDisp.resize(m_mesh[iMesh]->numPoints()*3,0.); |
125 |
|
✗ |
io[iMesh]->readVariable(1, iteration,&m_vectorDisp[0], m_vectorDisp.size()); |
126 |
|
|
} |
127 |
|
✗ |
std::transform(m_vectorDisp.begin(), m_vectorDisp.end(), m_vectorDisp.begin(),std::bind(std::multiplies<double>(),FelisceParam::instance().spaceUnit, std::placeholders::_1)); |
128 |
|
|
} |
129 |
|
|
|
130 |
|
|
|
131 |
|
✗ |
void LinearProblemBidomain::writeEnsightScalar(double* solValue, int idIter, std::string varName) { |
132 |
|
✗ |
std::string iteration; |
133 |
|
✗ |
if (idIter < 10) |
134 |
|
✗ |
iteration = "0000" + std::to_string(idIter); |
135 |
|
✗ |
else if (idIter < 100) |
136 |
|
✗ |
iteration = "000" + std::to_string(idIter); |
137 |
|
✗ |
else if (idIter < 1000) |
138 |
|
✗ |
iteration = "00" + std::to_string(idIter); |
139 |
|
✗ |
else if (idIter < 10000) |
140 |
|
✗ |
iteration = "0" + std::to_string(idIter); |
141 |
|
✗ |
else if (idIter < 100000) |
142 |
|
✗ |
iteration = std::to_string(idIter); |
143 |
|
|
|
144 |
|
✗ |
std::string fileName = FelisceParam::instance().resultDir + "/" + varName + "." + iteration + ".scl"; |
145 |
|
|
FILE * solFile; |
146 |
|
✗ |
solFile = fopen(fileName.c_str(),"w"); |
147 |
|
✗ |
fprintf(solFile, "Scalar per node\n"); |
148 |
|
✗ |
int count = 0; |
149 |
|
✗ |
for(int j=0; j < static_cast<int>(m_numDof/2); j++) { |
150 |
|
✗ |
fprintf(solFile,"%12.5e", solValue[j]); |
151 |
|
✗ |
count++; |
152 |
|
✗ |
if(count == 6) { |
153 |
|
✗ |
fprintf(solFile,"\n"); |
154 |
|
✗ |
count = 0; |
155 |
|
|
} |
156 |
|
|
} |
157 |
|
✗ |
fprintf(solFile,"\n"); |
158 |
|
✗ |
fclose(solFile); |
159 |
|
|
} |
160 |
|
|
|
161 |
|
✗ |
double LinearProblemBidomain::conductivityHom(int currentLabel) { |
162 |
|
✗ |
double conductivity = 1.0; |
163 |
|
|
|
164 |
|
✗ |
if (FelisceParam::instance().hasHeteroCondAtria) { |
165 |
|
✗ |
if (std::fabs(currentLabel - 31.0) < 1.0e-12) { //SN |
166 |
|
✗ |
conductivity = 0.3; |
167 |
|
✗ |
} else if (std::fabs(currentLabel - 32.0) < 1.0e-12) { //CT |
168 |
|
✗ |
conductivity = 3.0; |
169 |
|
✗ |
} else if (std::fabs(currentLabel - 33.0) < 1.0e-12) { //BB |
170 |
|
✗ |
conductivity = 4.75; |
171 |
|
✗ |
} else if (std::fabs(currentLabel - 34.0) < 1.0e-12) { //RAI |
172 |
|
✗ |
conductivity = 0.375; |
173 |
|
✗ |
} else if (std::fabs(currentLabel - 35.0) < 1.0e-12) { //PM |
174 |
|
✗ |
conductivity = 1.5; |
175 |
|
|
} else { //regular atria and F0 |
176 |
|
✗ |
conductivity = 1.0; |
177 |
|
|
} |
178 |
|
|
} |
179 |
|
|
|
180 |
|
✗ |
return conductivity; |
181 |
|
|
} |
182 |
|
|
|
183 |
|
|
|
184 |
|
✗ |
double LinearProblemBidomain::conductivityHet(int currentLabel) { |
185 |
|
✗ |
double conductivity = 1.0; |
186 |
|
|
|
187 |
|
✗ |
if (FelisceParam::instance().hasHeteroCondAtria) { |
188 |
|
✗ |
if (std::fabs(currentLabel - 31.0) < 1.0e-12) { //SN |
189 |
|
✗ |
conductivity = 0.36; |
190 |
|
✗ |
} else if (std::fabs(currentLabel - 32.0) < 1.0e-12) { //CT |
191 |
|
✗ |
conductivity = 4.5; |
192 |
|
✗ |
} else if (std::fabs(currentLabel - 33.0) < 1.0e-12) { //BB |
193 |
|
✗ |
conductivity = 7.75; |
194 |
|
✗ |
} else if (std::fabs(currentLabel - 34.0) < 1.0e-12) { //RAI |
195 |
|
✗ |
conductivity = 0.45; |
196 |
|
✗ |
} else if (std::fabs(currentLabel - 35.0) < 1.0e-12) { //PM |
197 |
|
✗ |
conductivity = 1.8; |
198 |
|
✗ |
} else if (std::fabs(currentLabel - 36.0) < 1.0e-12) { //F0 |
199 |
|
✗ |
conductivity = 0.9; |
200 |
|
|
} else { //regular atria |
201 |
|
✗ |
conductivity = 1.0; |
202 |
|
|
} |
203 |
|
|
} |
204 |
|
|
|
205 |
|
✗ |
return conductivity; |
206 |
|
|
} |
207 |
|
|
|
208 |
|
|
|
209 |
|
✗ |
void LinearProblemBidomain::writeEnsightCase(int numIt, std::string varName) { |
210 |
|
|
|
211 |
|
|
FILE * pFile; |
212 |
|
✗ |
std::string caseName = FelisceParam::instance().resultDir + "/" + varName + ".case"; |
213 |
|
✗ |
pFile = fopen(caseName.c_str(),"w"); |
214 |
|
✗ |
fprintf(pFile,"FORMAT\n"); |
215 |
|
✗ |
fprintf(pFile,"type: ensight\n"); |
216 |
|
✗ |
fprintf(pFile, "GEOMETRY\n"); |
217 |
|
✗ |
std::string nameGeo = "model: 1 " + FelisceParam::instance().outputMesh[0] + "\n"; |
218 |
|
✗ |
fprintf(pFile,"%s", nameGeo.c_str()); |
219 |
|
✗ |
fprintf(pFile, "VARIABLE\n"); |
220 |
|
✗ |
std::string secName = "scalar per node: 1 " + varName + " " + varName +".*****.scl\n"; |
221 |
|
✗ |
fprintf(pFile, "%s", secName.c_str()); |
222 |
|
✗ |
fprintf(pFile, "TIME\n"); |
223 |
|
✗ |
fprintf(pFile, "time std::set: %d \n", 1); |
224 |
|
✗ |
fprintf(pFile, "number of steps: %d \n", numIt); |
225 |
|
✗ |
fprintf(pFile, "filename start number: %d\n", 0); |
226 |
|
✗ |
fprintf(pFile, "filename increment: %d\n", 1); |
227 |
|
✗ |
fprintf(pFile, "time values:\n"); |
228 |
|
|
|
229 |
|
✗ |
const double dt = FelisceParam::instance().timeStep * FelisceParam::instance().frequencyWriteSolution; |
230 |
|
✗ |
int count=0; |
231 |
|
✗ |
for(int j=0; j < numIt; j++) { |
232 |
|
✗ |
double valT = j*dt; |
233 |
|
✗ |
fprintf(pFile,"%lf ",valT); |
234 |
|
✗ |
count++; |
235 |
|
✗ |
if(count == 6) { |
236 |
|
✗ |
fprintf(pFile,"\n"); |
237 |
|
✗ |
count=0; |
238 |
|
|
} |
239 |
|
|
} |
240 |
|
✗ |
fclose(pFile); |
241 |
|
|
} |
242 |
|
|
|
243 |
|
|
|
244 |
|
✗ |
void LinearProblemBidomain::writeEnsightCaseCurrents(int numIt, std::vector<std::string> varName) { |
245 |
|
|
|
246 |
|
|
FILE * pFile; |
247 |
|
✗ |
std::string caseName = FelisceParam::instance().resultDir + "/" + "courants" + ".case"; |
248 |
|
✗ |
pFile = fopen(caseName.c_str(),"w"); |
249 |
|
✗ |
fprintf(pFile,"FORMAT\n"); |
250 |
|
✗ |
fprintf(pFile,"type: ensight\n"); |
251 |
|
✗ |
fprintf(pFile, "GEOMETRY\n"); |
252 |
|
✗ |
std::string nameGeo = "model: 1 " + FelisceParam::instance().outputMesh[currentMesh()] + "\n"; |
253 |
|
✗ |
fprintf(pFile,"%s", nameGeo.c_str()); |
254 |
|
✗ |
fprintf(pFile, "VARIABLE\n"); |
255 |
|
|
|
256 |
|
✗ |
const double nbvarName = varName.size(); |
257 |
|
✗ |
for(int i=0; i<nbvarName; i++) { |
258 |
|
✗ |
std::string secName = "scalar per node: 1 " + varName[i] + " " + varName[i] +".*****.scl\n"; |
259 |
|
✗ |
fprintf(pFile, "%s", secName.c_str()); |
260 |
|
|
} |
261 |
|
|
|
262 |
|
✗ |
fprintf(pFile, "TIME\n"); |
263 |
|
✗ |
fprintf(pFile, "time std::set: %d \n", 1); |
264 |
|
✗ |
fprintf(pFile, "number of steps: %d \n", numIt); |
265 |
|
✗ |
fprintf(pFile, "filename start number: %d\n", 0); |
266 |
|
✗ |
fprintf(pFile, "filename increment: %d\n", 1); |
267 |
|
✗ |
fprintf(pFile, "time values:\n"); |
268 |
|
|
|
269 |
|
✗ |
const double dt = FelisceParam::instance().timeStep * FelisceParam::instance().frequencyWriteSolution; |
270 |
|
✗ |
int count=0; |
271 |
|
✗ |
for(int j=0; j < numIt; j++) { |
272 |
|
✗ |
double valT = j*dt; |
273 |
|
✗ |
fprintf(pFile,"%lf ",valT); |
274 |
|
✗ |
count++; |
275 |
|
✗ |
if(count == 6) { |
276 |
|
✗ |
fprintf(pFile,"\n"); |
277 |
|
✗ |
count=0; |
278 |
|
|
} |
279 |
|
|
} |
280 |
|
✗ |
fclose(pFile); |
281 |
|
|
} |
282 |
|
|
|
283 |
|
|
|
284 |
|
✗ |
void LinearProblemBidomain::writeSolution(int rank, std::vector<IO::Pointer>& io, double& time, int iteration) { |
285 |
|
✗ |
LinearProblem::writeSolution(rank,io, time, iteration); |
286 |
|
|
// Add fiber as an output variable. |
287 |
|
|
/* |
288 |
|
|
io.writeSolution(rank, time, iteration, 1, "fiber", m_fiber, m_supportDofUnknown[0].listNode().size(), m_dimension); |
289 |
|
|
if (FelisceParam::instance().typeOfAppliedCurrent == "zygote") |
290 |
|
|
{ |
291 |
|
|
// Add distance std::unordered_map as an output variable. |
292 |
|
|
io.writeSolution(rank, time, iteration, 0, "distance", m_endocardiumDistance, m_supportDofUnknown[0].listNode().size(), m_dimension); |
293 |
|
|
} |
294 |
|
|
*/ |
295 |
|
|
} |
296 |
|
|
|
297 |
|
|
|
298 |
|
✗ |
void LinearProblemBidomain::getFiberDirection(felInt iel, int iUnknown, std::vector<double>& elemFiber) { |
299 |
|
|
//Warning: the fibers must be normalized. |
300 |
|
✗ |
int numSupport = static_cast<int>(supportDofUnknown(iUnknown).iEle()[iel+1] - supportDofUnknown(iUnknown).iEle()[iel]); |
301 |
|
✗ |
elemFiber.resize( numSupport*3,0.); |
302 |
|
✗ |
for (int i = 0; i < numSupport; i++) { |
303 |
|
✗ |
elemFiber[i*3] = m_vectorFiber[3*(supportDofUnknown(iUnknown).iSupportDof()[supportDofUnknown(iUnknown).iEle()[iel]+i])]; |
304 |
|
✗ |
elemFiber[i*3+1] = m_vectorFiber[3*(supportDofUnknown(iUnknown).iSupportDof()[supportDofUnknown(iUnknown).iEle()[iel]+i])+1]; |
305 |
|
✗ |
elemFiber[i*3+2] = m_vectorFiber[3*(supportDofUnknown(iUnknown).iSupportDof()[supportDofUnknown(iUnknown).iEle()[iel]+i])+2]; |
306 |
|
|
} |
307 |
|
|
} |
308 |
|
|
|
309 |
|
✗ |
void LinearProblemBidomain::getAngle(felInt iel, int iUnknown, std::vector<double>& elemAngle) { |
310 |
|
✗ |
int numSupport = static_cast<int>(supportDofUnknown(iUnknown).iEle()[iel+1] - supportDofUnknown(iUnknown).iEle()[iel]); |
311 |
|
✗ |
elemAngle.resize(numSupport,0.); |
312 |
|
✗ |
for (int i = 0; i < numSupport; i++) { |
313 |
|
✗ |
elemAngle[i] = m_vectorAngle[(supportDofUnknown(iUnknown).iSupportDof()[supportDofUnknown(iUnknown).iEle()[iel]+i])]; |
314 |
|
|
} |
315 |
|
|
} |
316 |
|
|
|
317 |
|
✗ |
void LinearProblemBidomain::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS) { |
318 |
|
|
IGNORE_UNUSED_ELT_TYPE; |
319 |
|
|
IGNORE_UNUSED_FLAG_MATRIX_RHS; |
320 |
|
|
//Init pointer on Finite Element, Variable or idVariable |
321 |
|
✗ |
m_ipotTransMemb = m_listVariable.getVariableIdList(potTransMemb); |
322 |
|
✗ |
m_ipotExtraCell = m_listVariable.getVariableIdList(potExtraCell); |
323 |
|
✗ |
m_fePotTransMemb = m_listCurrentFiniteElement[m_ipotTransMemb]; |
324 |
|
✗ |
m_fePotExtraCell = m_fePotTransMemb; |
325 |
|
|
|
326 |
|
|
} |
327 |
|
|
|
328 |
|
✗ |
void LinearProblemBidomain::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) { |
329 |
|
|
IGNORE_UNUSED_FLAG_MATRIX_RHS; |
330 |
|
|
IGNORE_UNUSED_ELEM_ID_POINT; |
331 |
|
✗ |
m_fePotTransMemb->updateMeasQuadPt(0, elemPoint); |
332 |
|
✗ |
m_fePotTransMemb->updateFirstDeriv(0, elemPoint); |
333 |
|
|
|
334 |
|
|
// Assembling matrices[0] |
335 |
|
✗ |
double coefTime = 1./m_fstransient->timeStep; |
336 |
|
✗ |
double coefAm = FelisceParam::instance().Am; |
337 |
|
✗ |
double coefCm = FelisceParam::instance().Cm; |
338 |
|
✗ |
double coef = m_bdf->coeffDeriv0()*coefTime*coefAm*coefCm; |
339 |
|
|
|
340 |
|
|
//1st equation : Am*Cm/dt V_m |
341 |
|
✗ |
m_elementMat[0]->phi_i_phi_j(coef,*m_fePotTransMemb,0,0,1); |
342 |
|
|
|
343 |
|
✗ |
if (FelisceParam::instance().testCase == 1) { |
344 |
|
✗ |
std::vector<double> elemFiber; |
345 |
|
✗ |
getFiberDirection(iel,m_ipotTransMemb, elemFiber); |
346 |
|
|
|
347 |
|
✗ |
if (FelisceParam::instance().model == "Bidomain") { |
348 |
|
|
//1st equation : \sigma_i^t \grad V_m |
349 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1); |
350 |
|
|
|
351 |
|
|
//1st equation : (\sigma_i^l-\sigma_i^t) a vec a \grad V_m |
352 |
|
✗ |
m_elementMat[0]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1); |
353 |
|
|
|
354 |
|
|
//1st equation : \sigma_i^t \grad u_e |
355 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotExtraCell,0,1,1); |
356 |
|
|
|
357 |
|
|
//1st equation : (\sigma_i^l-\sigma_i^t) a vec a \grad u_e |
358 |
|
✗ |
m_elementMat[0]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotExtraCell,0,1,1); |
359 |
|
|
|
360 |
|
|
//2nd equation : \sigma_i^t \grad V_m |
361 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,1,0,1); |
362 |
|
|
|
363 |
|
|
//2nd equation : (\sigma_i^l-\sigma_i^t) a vec a \grad V_m |
364 |
|
✗ |
m_elementMat[0]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,1,0,1); |
365 |
|
|
|
366 |
|
|
//2nd equation : (\sigma_i^t+\sigma_e^t) \grad u_e |
367 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraTransvTensor,*m_fePotExtraCell,1,1,1); |
368 |
|
|
|
369 |
|
|
//2nd equation : ((\sigma_i^l-\sigma_i^t)+(\sigma_e^l-\sigma_e^t)) a vec a \grad u_e |
370 |
|
✗ |
m_elementMat[0]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraFiberTensor-FelisceParam::instance().extraTransvTensor,elemFiber,*m_fePotExtraCell,1,1,1); |
371 |
|
|
|
372 |
|
|
} |
373 |
|
✗ |
else if (FelisceParam::instance().model == "Monodomain") { |
374 |
|
|
//1st equation : \sigma_i^t \grad V_m |
375 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1); |
376 |
|
|
|
377 |
|
|
//1st equation : (\sigma_i^l-\sigma_i^t) a vec a \grad V_m |
378 |
|
✗ |
m_elementMat[0]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1); |
379 |
|
|
|
380 |
|
|
//2nd equation : (\sigma_i^t+\sigma_e^t) \grad u_e |
381 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraTransvTensor,*m_fePotExtraCell,1,1,1); |
382 |
|
|
|
383 |
|
|
//2nd equation : ((\sigma_i^l-\sigma_i^t)+(\sigma_e^l-\sigma_e^t)) a vec a \grad u_e |
384 |
|
✗ |
m_elementMat[0]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraFiberTensor-FelisceParam::instance().extraTransvTensor,elemFiber,*m_fePotExtraCell,1,1,1); |
385 |
|
|
} |
386 |
|
|
} |
387 |
|
|
else { |
388 |
|
✗ |
if (FelisceParam::instance().model == "Bidomain") { |
389 |
|
|
//1st equation : \sigma_i^t \grad V_m |
390 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1); |
391 |
|
|
|
392 |
|
|
//1st equation : \sigma_i^t \grad u_e |
393 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotExtraCell,0,1,1); |
394 |
|
|
|
395 |
|
|
//2nd equation : \sigma_i^t \grad V_m |
396 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,1,0,1); |
397 |
|
|
|
398 |
|
|
//2nd equation : (\sigma_i^t+\sigma_e^t) \grad u_e |
399 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraTransvTensor,*m_fePotExtraCell,1,1,1); |
400 |
|
|
|
401 |
|
|
} |
402 |
|
✗ |
else if (FelisceParam::instance().model == "Monodomain") { |
403 |
|
|
//1st equation : \sigma_i^t \grad V_m |
404 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1); |
405 |
|
|
|
406 |
|
|
//2nd equation : (\sigma_i^t+\sigma_e^t) \grad u_e |
407 |
|
✗ |
m_elementMat[0]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraTransvTensor,*m_fePotExtraCell,1,1,1); |
408 |
|
|
} |
409 |
|
|
} |
410 |
|
|
//Assembling matrix(1) = mass (only first quadrant). |
411 |
|
✗ |
m_elementMat[1]->phi_i_phi_j(1.,*m_fePotTransMemb,0,0,1); |
412 |
|
✗ |
m_elementMat[1]->phi_i_phi_j(1.,*m_fePotTransMemb,1,1,1); |
413 |
|
|
|
414 |
|
✗ |
if (FelisceParam::instance().stateFilter) { |
415 |
|
✗ |
*m_elementMat[2]=*m_elementMat[0]; |
416 |
|
|
} |
417 |
|
|
} |
418 |
|
|
|
419 |
|
✗ |
void LinearProblemBidomain::computeElementArraySurfaceModel(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) { |
420 |
|
|
(void) elemIdPoint; |
421 |
|
|
(void) flagMatrixRHS; |
422 |
|
✗ |
m_ipotTransMemb = m_listVariable.getVariableIdList(potTransMemb); |
423 |
|
✗ |
m_fePotTransMembCurv = m_listCurvilinearFiniteElement[m_ipotTransMemb]; |
424 |
|
✗ |
m_fePotExtraCellCurv = m_fePotTransMembCurv; |
425 |
|
|
|
426 |
|
✗ |
m_fePotTransMembCurv->updateMeasNormalContra(0, elemPoint); |
427 |
|
|
|
428 |
|
✗ |
auto bndRefIt = find(FelisceParam::instance().notBoundarySurfaceLabels.begin(), FelisceParam::instance().notBoundarySurfaceLabels.end(), m_currentLabel); |
429 |
|
|
|
430 |
|
✗ |
if(bndRefIt != FelisceParam::instance().notBoundarySurfaceLabels.end()) { |
431 |
|
|
// if ((m_currentLabel > 29) && (m_currentLabel < 41)) { |
432 |
|
✗ |
std::vector<double> elemAngle; |
433 |
|
✗ |
getAngle(iel,m_ipotTransMemb,elemAngle); |
434 |
|
|
|
435 |
|
✗ |
std::vector<double> elemFiber; |
436 |
|
✗ |
getFiberDirection(iel,m_ipotTransMemb,elemFiber); |
437 |
|
|
|
438 |
|
|
// Assembling matrix(0) |
439 |
|
✗ |
double coefTime = 1./m_fstransient->timeStep; |
440 |
|
✗ |
double coefAm = FelisceParam::instance().Am; |
441 |
|
✗ |
double coefCm = FelisceParam::instance().Cm; |
442 |
|
✗ |
double coef = m_bdf->coeffDeriv0()*coefTime*coefAm*coefCm; |
443 |
|
|
|
444 |
|
✗ |
double atriaThick=0.2; // cm - atria thickness |
445 |
|
|
|
446 |
|
|
//1st equation : Am*Cm/dt V_m |
447 |
|
✗ |
m_elementMatBD[0]->phi_i_phi_j(atriaThick*coef,*m_fePotTransMembCurv,0,0,1); |
448 |
|
|
|
449 |
|
|
//1st equation : \sigma_i^t \grad V_m |
450 |
|
✗ |
if (m_currentLabel < 40) { // regular atrial myocardium |
451 |
|
✗ |
m_elementMatBD[0]->grad_phi_i_grad_phi_j(atriaThick*FelisceParam::instance().intraTransvTensorAtria*conductivityHom(m_currentLabel),*m_fePotTransMembCurv,0,0,1); |
452 |
|
✗ |
m_elementMatBD[0]->tau_orthotau_grad_phi_i_grad_phi_j(atriaThick*(FelisceParam::instance().intraFiberTensorAtria - FelisceParam::instance().intraTransvTensorAtria)*conductivityHet(m_currentLabel),elemFiber,elemAngle,*m_fePotTransMembCurv,0,0,1); |
453 |
|
|
} |
454 |
|
|
|
455 |
|
|
//1st equation : \sigma_i^t \grad u_e |
456 |
|
✗ |
if (m_currentLabel < 40) { // regular atrial myocardium |
457 |
|
✗ |
m_elementMatBD[0]->grad_phi_i_grad_phi_j(atriaThick*FelisceParam::instance().intraTransvTensorAtria*conductivityHom(m_currentLabel),*m_fePotExtraCellCurv,0,1,1); |
458 |
|
✗ |
m_elementMatBD[0]->tau_orthotau_grad_phi_i_grad_phi_j(atriaThick*(FelisceParam::instance().intraFiberTensorAtria - FelisceParam::instance().intraTransvTensorAtria)*conductivityHet(m_currentLabel),elemFiber,elemAngle,*m_fePotExtraCellCurv,0,1,1); |
459 |
|
|
} |
460 |
|
|
|
461 |
|
|
//2nd equation : \sigma_i^t \grad V_m |
462 |
|
✗ |
if (m_currentLabel < 40) { // regular atrial myocardium |
463 |
|
✗ |
m_elementMatBD[0]->grad_phi_i_grad_phi_j(atriaThick*FelisceParam::instance().intraTransvTensorAtria*conductivityHom(m_currentLabel),*m_fePotTransMembCurv,1,0,1); |
464 |
|
✗ |
m_elementMatBD[0]->tau_orthotau_grad_phi_i_grad_phi_j(atriaThick*(FelisceParam::instance().intraFiberTensorAtria - FelisceParam::instance().intraTransvTensorAtria)*conductivityHet(m_currentLabel),elemFiber,elemAngle,*m_fePotTransMembCurv,1,0,1); |
465 |
|
|
} |
466 |
|
|
|
467 |
|
|
//2nd equation : (\sigma_i^t+\sigma_e^t) \grad u_e |
468 |
|
✗ |
double connectionCoeff = 3.e-4; // 0.0003 |
469 |
|
✗ |
if (m_currentLabel < 40) { // atrial myocardium |
470 |
|
✗ |
m_elementMatBD[0]->grad_phi_i_grad_phi_j(atriaThick*(FelisceParam::instance().intraTransvTensorAtria+FelisceParam::instance().extraTransvTensorAtria)*conductivityHom(m_currentLabel),*m_fePotExtraCellCurv,1,1,1); |
471 |
|
✗ |
m_elementMatBD[0]->tau_orthotau_grad_phi_i_grad_phi_j(atriaThick*(FelisceParam::instance().intraFiberTensorAtria - FelisceParam::instance().intraTransvTensorAtria + FelisceParam::instance().extraFiberTensorAtria - FelisceParam::instance().extraTransvTensorAtria)*conductivityHet(m_currentLabel),elemFiber,elemAngle,*m_fePotExtraCellCurv,1,1,1); |
472 |
|
✗ |
} else if (std::fabs(m_currentLabel - 40) < 1.0e-8) { // connection aera |
473 |
|
✗ |
m_elementMatBD[0]->grad_phi_i_grad_phi_j(atriaThick*(FelisceParam::instance().extraTransvTensorAtria)*connectionCoeff,*m_fePotExtraCellCurv,1,1,1); |
474 |
|
|
} |
475 |
|
|
|
476 |
|
|
//Assembling matrix(1) = mass (only first quadrant if hasAppliedExteriorCurrent = false) |
477 |
|
✗ |
m_elementMatBD[1]->phi_i_phi_j(atriaThick*1.,*m_fePotTransMembCurv,0,0,1); |
478 |
|
✗ |
m_elementMatBD[1]->phi_i_phi_j(atriaThick*1.,*m_fePotTransMembCurv,1,1,1); |
479 |
|
|
} |
480 |
|
|
} |
481 |
|
|
|
482 |
|
✗ |
void LinearProblemBidomain::computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, felInt& iel, FlagMatrixRHS /*flagMatrixRHS*/) { |
483 |
|
|
|
484 |
|
|
(void)iel; |
485 |
|
|
|
486 |
|
✗ |
m_curvFe->updateMeas(0,elemPoint); |
487 |
|
|
// #ifdef FELISCE_WITH_CVGRAPH |
488 |
|
|
// if (FelisceParam::instance().withCVG ) { |
489 |
|
|
// this->cvgraphNaturalBC(iel); |
490 |
|
|
// } |
491 |
|
|
// #endif |
492 |
|
|
} |
493 |
|
|
|
494 |
|
✗ |
void LinearProblemBidomain::addMatrixRHS() { |
495 |
|
|
|
496 |
|
|
//Rebuild the matrix with the filter |
497 |
|
✗ |
if (FelisceParam::instance().stateFilter ) { |
498 |
|
✗ |
matrix(0).copyFrom(matrix(2),SAME_NONZERO_PATTERN); |
499 |
|
|
} |
500 |
|
|
|
501 |
|
|
// _RHS = m_mass * _RHS. |
502 |
|
✗ |
PetscVector tmpB; |
503 |
|
✗ |
tmpB.duplicateFrom(vector()); |
504 |
|
✗ |
tmpB.copyFrom(vector()); |
505 |
|
✗ |
mult(matrix(1),tmpB,vector()); // _RHS = m_mass * tmpB |
506 |
|
✗ |
tmpB.destroy(); |
507 |
|
|
|
508 |
|
✗ |
if (m_initializeCVGraphInterface) { |
509 |
|
✗ |
vector(1).copyFrom(vector()); |
510 |
|
|
} |
511 |
|
|
|
512 |
|
✗ |
if(FelisceParam::instance().useROM) { |
513 |
|
✗ |
if (m_fstransient->iteration == 1) { |
514 |
|
✗ |
m_rom->fullToReducedMatrix(matrix(0)); |
515 |
|
|
} |
516 |
|
✗ |
m_rom->fullToReducedRHS(vector()); |
517 |
|
|
} |
518 |
|
|
|
519 |
|
|
} |
520 |
|
|
|
521 |
|
✗ |
void LinearProblemBidomain::readMatch() { |
522 |
|
|
// Input files. |
523 |
|
✗ |
std::string matchFileName = FelisceParam::instance().inputDirectory + "/" + FelisceParam::instance().ECGmatchFile;// load the match file (heart nodes) |
524 |
|
|
|
525 |
|
✗ |
std::ifstream matchFile(matchFileName.c_str()); |
526 |
|
✗ |
if(!matchFile) { |
527 |
|
✗ |
FEL_ERROR("Fatal error: cannot open match file " + matchFileName ); |
528 |
|
|
} |
529 |
|
✗ |
if (FelisceParam::verbose() > 0) { |
530 |
|
✗ |
std::cout << "Reading " << matchFileName << std::endl; |
531 |
|
|
} |
532 |
|
✗ |
matchFile >> m_numElectrodeNode; |
533 |
|
✗ |
m_electrodeNode = new felInt[m_numElectrodeNode]; |
534 |
|
|
felInt value; |
535 |
|
✗ |
for(int j=0; j<m_numElectrodeNode; j++) { |
536 |
|
✗ |
matchFile >> value; |
537 |
|
✗ |
m_electrodeNode[j] = value; |
538 |
|
|
} |
539 |
|
✗ |
matchFile.close(); |
540 |
|
|
} |
541 |
|
|
|
542 |
|
|
// Elisa, Dec. 2014 - read measures for Luenberger filter |
543 |
|
✗ |
void LinearProblemBidomain::readElectrodeMeasure() { |
544 |
|
|
// load the measures |
545 |
|
✗ |
std::string inputFileName = FelisceParam::instance().inputDirectory + "/" + FelisceParam::instance().ECGfileName + ".tab"; |
546 |
|
✗ |
std::string outputFileName = FelisceParam::instance().resultDir + "/" + FelisceParam::instance().ECGfileName + "_control.tab"; |
547 |
|
✗ |
int ierr = 0; |
548 |
|
✗ |
if ( !filesystemUtil::directoryExists(FelisceParam::instance().resultDir.c_str()) ) { |
549 |
|
✗ |
std::string command = "mkdir -p " + FelisceParam::instance().resultDir; |
550 |
|
✗ |
ierr = system(command.c_str()); |
551 |
|
✗ |
if (ierr) { |
552 |
|
✗ |
std::cout << "Cannot copy control file in resultDir." << std::endl; |
553 |
|
|
} |
554 |
|
|
} |
555 |
|
✗ |
if (ierr == 0) { |
556 |
|
✗ |
filesystemUtil::copyFileSomewhereElse(inputFileName,outputFileName); |
557 |
|
|
} |
558 |
|
|
|
559 |
|
✗ |
std::ifstream ecgFile(inputFileName.c_str(), std::ios_base::in); |
560 |
|
✗ |
if ( !ecgFile.is_open() ) { |
561 |
|
✗ |
FEL_ERROR(" ERROR: Can not open file " + inputFileName); |
562 |
|
|
} |
563 |
|
|
|
564 |
|
✗ |
m_electrodeMeasure = new double* [m_numElectrodeNode+1]; |
565 |
|
✗ |
for (int i=0; i<m_numElectrodeNode+1; i++) { |
566 |
|
✗ |
m_electrodeMeasure[i] = new double [FelisceParam::instance().timeIterationMax]; |
567 |
|
|
} |
568 |
|
|
|
569 |
|
✗ |
for (int j=0; j<FelisceParam::instance().timeIterationMax; j++) { |
570 |
|
✗ |
for(int i=0; i<m_numElectrodeNode+1; i++) { |
571 |
|
✗ |
ecgFile >> m_electrodeMeasure[i][j]; |
572 |
|
|
} |
573 |
|
|
} |
574 |
|
✗ |
ecgFile.close(); |
575 |
|
|
|
576 |
|
✗ |
m_electrodeControl.duplicateFrom(vector()); |
577 |
|
|
} |
578 |
|
|
|
579 |
|
|
// Elisa, Dec. 2014 - Luenberger filter |
580 |
|
✗ |
void LinearProblemBidomain::addElectrodeCondtrol() { |
581 |
|
✗ |
if (m_electrodeNode == nullptr) { |
582 |
|
✗ |
readMatch(); |
583 |
|
|
} |
584 |
|
✗ |
if (m_electrodeMeasure == nullptr) { |
585 |
|
✗ |
readElectrodeMeasure(); |
586 |
|
|
} |
587 |
|
|
|
588 |
|
✗ |
if (FelisceParam::instance().electrodeControl) { |
589 |
|
✗ |
if ( !Tools::almostEqual(m_electrodeMeasure[0][m_fstransient->iteration-1], m_fstransient->time-m_fstransient->timeStep, m_fstransient->timeStep/10.) ) { |
590 |
|
✗ |
FEL_WARNING("Measures physical time different from solver physical time.\n"); |
591 |
|
|
} |
592 |
|
|
} |
593 |
|
|
|
594 |
|
✗ |
m_electrodeControl.set(0.); |
595 |
|
|
|
596 |
|
✗ |
int sizePot = numDofLocalPerUnknown(potExtraCell); |
597 |
|
✗ |
double value[sizePot]; |
598 |
|
✗ |
felInt idLocalUe[sizePot]; |
599 |
|
✗ |
felInt idGlobalUe[sizePot]; |
600 |
|
✗ |
felInt idLocalVm[sizePot]; |
601 |
|
✗ |
felInt idGlobalVm[sizePot]; |
602 |
|
|
|
603 |
|
✗ |
for (felInt i = 0; i < sizePot; i++) { |
604 |
|
✗ |
idLocalUe[i] = i; |
605 |
|
✗ |
idLocalVm[i] = i; |
606 |
|
✗ |
value[i] = 0.; |
607 |
|
|
} |
608 |
|
✗ |
ISLocalToGlobalMappingApply(mappingDofLocalToDofGlobal(potTransMemb),sizePot,&idLocalVm[0],&idGlobalVm[0]); |
609 |
|
✗ |
ISLocalToGlobalMappingApply(mappingDofLocalToDofGlobal(potExtraCell),sizePot,&idLocalUe[0],&idGlobalUe[0]); |
610 |
|
|
|
611 |
|
✗ |
for (felInt i = 0; i < m_numElectrodeNode; i++) { |
612 |
|
✗ |
felInt idElectrode = m_electrodeNode[i]-1; |
613 |
|
✗ |
value[idGlobalVm[idElectrode]] = m_electrodeMeasure[i+1][m_fstransient->iteration-1]; |
614 |
|
|
double valuePot; |
615 |
|
✗ |
this->solution().getValues(1,&idGlobalUe[idElectrode],&valuePot); |
616 |
|
✗ |
value[idGlobalVm[idElectrode]] -= valuePot; |
617 |
|
|
} |
618 |
|
|
|
619 |
|
✗ |
m_electrodeControl.setValues(sizePot,&idGlobalVm[0],&value[0],INSERT_VALUES); |
620 |
|
|
|
621 |
|
✗ |
m_electrodeControl.assembly(); |
622 |
|
|
|
623 |
|
✗ |
vector().axpy(-FelisceParam::instance().aFilter,m_electrodeControl); |
624 |
|
|
|
625 |
|
✗ |
vector().assembly(); |
626 |
|
|
|
627 |
|
|
} |
628 |
|
|
|
629 |
|
✗ |
void LinearProblemBidomain::assembleOnlyMassMatrix(PetscMatrix& mass) |
630 |
|
|
{ |
631 |
|
|
// Allocation |
632 |
|
✗ |
mass.duplicateFrom(matrix(0),MAT_DO_NOT_COPY_VALUES); |
633 |
|
✗ |
mass.zeroEntries(); |
634 |
|
|
|
635 |
|
|
/* Loop function */ |
636 |
|
✗ |
auto function = [this,&mass](ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, felInt& ielSupportDof) { |
637 |
|
✗ |
setElemPoint(eltType, iel, elemPoint, elemIdPoint, &ielSupportDof); |
638 |
|
✗ |
m_elementMat[0]->zero(); |
639 |
|
✗ |
userComputeMass(elemPoint, elemIdPoint, ielSupportDof); |
640 |
|
✗ |
setValueCustomMatrix(ielSupportDof, ielSupportDof, mass); |
641 |
|
✗ |
return 0.0; |
642 |
|
✗ |
}; |
643 |
|
|
/* Pre-function */ |
644 |
|
✗ |
auto functionpre = [this](ElementType& eltType) { |
645 |
|
✗ |
defineFiniteElement(eltType); |
646 |
|
✗ |
initElementArray(); |
647 |
|
✗ |
allocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix); |
648 |
|
✗ |
}; |
649 |
|
|
/* Post-function */ |
650 |
|
✗ |
auto functionpost = [this](ElementType& eltType) { |
651 |
|
|
IGNORE_UNUSED_ARGUMENT(eltType); |
652 |
|
✗ |
desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix); |
653 |
|
✗ |
}; |
654 |
|
✗ |
FEUtilities::LoopOverElements(*m_meshLocal[m_currentMesh], m_meshLocal[m_currentMesh]->bagElementTypeDomain(), &function, &functionpre, &functionpost); |
655 |
|
|
|
656 |
|
|
// Final assembly |
657 |
|
✗ |
mass.assembly(MAT_FINAL_ASSEMBLY); |
658 |
|
|
} |
659 |
|
|
|
660 |
|
✗ |
void LinearProblemBidomain::initAndStoreMassMatrix() { |
661 |
|
✗ |
assembleOnlyMassMatrix(m_massMatrix); |
662 |
|
|
// TO DO: This should not be here, it's only temporary. Should be moved to prepareCVGRobin(). |
663 |
|
✗ |
initMassBoundaryForStimulation(); |
664 |
|
|
#ifdef FELISCE_WITH_CVGRAPH |
665 |
|
✗ |
if ( FelisceParam::instance().withCVG ) { |
666 |
|
✗ |
initMassBoundaryForCVG(); |
667 |
|
|
} |
668 |
|
|
#endif |
669 |
|
|
} |
670 |
|
|
|
671 |
|
|
#ifdef FELISCE_WITH_CVGRAPH |
672 |
|
✗ |
void LinearProblemBidomain::initMassBoundaryForCVG() { |
673 |
|
✗ |
this->assembleMassBoundaryOnly(&LinearProblemBidomain::massMatrixComputer, |
674 |
|
✗ |
this->slave()->interfaceLabels(/*iConn*/0), |
675 |
|
|
&LinearProblemBidomain::initPerETMass, |
676 |
|
|
&LinearProblemBidomain::updateFE, |
677 |
|
✗ |
this->dofBD(/*iConn*/0), this->massBDVec()[/*iConn*/0]); |
678 |
|
|
} |
679 |
|
|
|
680 |
|
|
/*! \brief Function to assemble the mass matrix |
681 |
|
|
* Function to be called in the assemblyLoopBoundaryGeneral |
682 |
|
|
*/ |
683 |
|
✗ |
void LinearProblemBidomain::massMatrixComputer(felInt ielSupportDof) { |
684 |
|
✗ |
this->m_elementMatBD[0]->zero(); |
685 |
|
✗ |
this->m_elementMatBD[0]->phi_i_phi_j(/*coef*/1.,*m_curvFe,1,1,1); |
686 |
|
✗ |
this->setValueMatrixBD(ielSupportDof); |
687 |
|
|
} |
688 |
|
|
void |
689 |
|
✗ |
LinearProblemBidomain::initPerETMass() { |
690 |
|
✗ |
felInt numDofTotal = 0; |
691 |
|
|
//pay attention this numDofTotal is bigger than expected since it counts for all the VARIABLES |
692 |
|
✗ |
for (std::size_t iFe = 0; iFe < this->m_listCurvilinearFiniteElement.size(); iFe++) {//this loop is on variables while it should be on unknown |
693 |
|
✗ |
numDofTotal += this->m_listCurvilinearFiniteElement[iFe]->numDof()*this->m_listVariable[iFe].numComponent(); |
694 |
|
|
} |
695 |
|
✗ |
m_globPosColumn.resize(numDofTotal); m_globPosRow.resize(numDofTotal); m_matrixValues.resize(numDofTotal*numDofTotal); |
696 |
|
|
|
697 |
|
✗ |
m_curvFe = this->m_listCurvilinearFiniteElement[ m_ipotExtraCell ]; |
698 |
|
|
} |
699 |
|
|
void |
700 |
|
✗ |
LinearProblemBidomain::updateFE(const std::vector<Point*>& elemPoint, const std::vector<int>&) { |
701 |
|
✗ |
m_curvFe->updateMeasNormal(0, elemPoint); |
702 |
|
|
} |
703 |
|
|
#endif |
704 |
|
|
|
705 |
|
|
#ifdef FELISCE_WITH_CVGRAPH |
706 |
|
✗ |
void LinearProblemBidomain::readData() { |
707 |
|
✗ |
for ( std::size_t iConn(0); iConn<this->slave()->numConnections(); ++iConn ) { |
708 |
|
|
//std::cout << "####### LinearProblemBidomain::readData(): " << this->slave()->numVarToRead(iConn) << " variables to read" << std::endl; |
709 |
|
✗ |
std::vector<PetscVector> vecs; |
710 |
|
✗ |
for ( std::size_t cVar(0); cVar<this->slave()->numVarToRead(iConn); ++cVar) { |
711 |
|
✗ |
std::string varToRead = this->slave()->readVariable(iConn,cVar); |
712 |
|
✗ |
if ( varToRead == "CURRENT" ) { |
713 |
|
✗ |
vecs.push_back(m_seqVecs.Get("cvgraphCURRENT")); |
714 |
|
✗ |
} else if ( varToRead == "POTENTIAL" ) { |
715 |
|
✗ |
vecs.push_back(m_seqVecs.Get("cvgraphPOTENTIAL")); |
716 |
|
|
} else { |
717 |
|
✗ |
std::cout<<"Connection: "<<iConn<<"--variable to read: "<< varToRead << std::endl; |
718 |
|
✗ |
FEL_ERROR("Error in variable to read"); |
719 |
|
|
} |
720 |
|
|
} |
721 |
|
✗ |
this->slave()->receiveData(vecs,iConn); |
722 |
|
|
} |
723 |
|
|
} |
724 |
|
✗ |
void LinearProblemBidomain::sendData() { |
725 |
|
|
//std::cout << "LinearProblemBidomain::sendData()" <<std::endl; |
726 |
|
✗ |
this->gatherSolution(); |
727 |
|
✗ |
for ( std::size_t iConn(0); iConn<this->slave()->numConnections(); ++iConn ) { |
728 |
|
✗ |
std::vector<PetscVector> vecs; |
729 |
|
✗ |
for ( std::size_t cVar(0); cVar<this->slave()->numVarToSend(iConn); ++cVar) { |
730 |
|
✗ |
std::string varToSend = this->slave()->sendVariable(iConn, cVar); |
731 |
|
✗ |
if ( varToSend == "POTENTIAL" ) { |
732 |
|
✗ |
vecs.push_back(this->sequentialSolution()); |
733 |
|
✗ |
} else if ( varToSend == "CURRENT" ) { |
734 |
|
|
// TODO: Add prepareResidual() and move this to startIterationCVG |
735 |
|
✗ |
this->computeResidual(); |
736 |
|
✗ |
vecs.push_back(this->seqResidual()); |
737 |
|
|
} else { |
738 |
|
|
//std::cout<<"Connection: "<<iConn<<"--variable to send: "<< varToSend << std::endl; |
739 |
|
✗ |
FEL_ERROR("Error in variable to read"); |
740 |
|
|
} |
741 |
|
|
} |
742 |
|
✗ |
this->slave()->sendData(vecs,iConn); |
743 |
|
|
} |
744 |
|
|
} |
745 |
|
✗ |
void LinearProblemBidomain::addCurrentToRHS() { |
746 |
|
|
//! for each dof on the boundary on this proc. |
747 |
|
✗ |
felInt numLocalDofInterface=this->dofBD(/*iBD*/0).numLocalDofInterface(); |
748 |
|
✗ |
std::vector<double> tmp(numLocalDofInterface); |
749 |
|
✗ |
m_seqVecs.Get("cvgraphCURRENT").getValues( numLocalDofInterface, this->dofBD(/*iBD*/0).loc2PetscVolPtr(), tmp.data() ); |
750 |
|
✗ |
this->vector().setValues(numLocalDofInterface,this->dofBD(/*iBD*/0).loc2PetscVolPtr(),tmp.data(),ADD_VALUES); |
751 |
|
✗ |
this->vector().assembly(); |
752 |
|
|
} |
753 |
|
|
|
754 |
|
✗ |
void LinearProblemBidomain::cvgAddRobinToRHS(){ |
755 |
|
|
//std::cout << "LinearProblemBidomain::cvgAddRobinToRHS()" << std::endl; |
756 |
|
|
//std::cout << " alpha_robin = " << FelisceParam::instance().alphaRobin[0/*iRobin*/] << std::endl; |
757 |
|
|
} |
758 |
|
|
|
759 |
|
|
#endif |
760 |
|
|
|
761 |
|
|
} |
762 |
|
|
|