Nektar++
FilterElectrogram.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: FilterElectrogram.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Outputs virtual electrograms at specific locations.
32//
33///////////////////////////////////////////////////////////////////////////////
34
37#include <iomanip>
38
39using namespace std;
40
41namespace Nektar
42{
45 "Electrogram", FilterElectrogram::create);
46
47/**
48 *
49 */
52 const std::shared_ptr<SolverUtils::EquationSystem> &pEquation,
53 const ParamMap &pParams)
54 : Filter(pSession, pEquation)
55{
56 // OutputFile
57 std::string ext = ".ecg";
58 m_outputFile = Filter::SetupOutput(ext, pParams);
59
60 // OutputFrequency
61 auto it = pParams.find("OutputFrequency");
62 if (it == pParams.end())
63 {
65 }
66 else
67 {
68 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
69 m_outputFrequency = floor(equ.Evaluate());
70 }
71
72 // Points
73 it = pParams.find("Points");
74 ASSERTL0(it != pParams.end(), "Missing parameter 'Points'.");
75 m_electrogramStream.str(it->second);
76 m_index = 0;
77}
78
79/**
80 *
81 */
83{
84}
85
86/**
87 *
88 */
91 const NekDouble &time)
92{
93 ASSERTL0(!m_electrogramStream.fail(), "No history points in stream.");
94
95 m_index = 0;
96 Array<OneD, NekDouble> gloCoord(3, 0.0);
97 LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
98
99 // Read electrogram points
100 // Always use dim = 3 to allow electrode to be above surface
101 const size_t dim = 3;
102 size_t i = 0;
103
104 while (!m_electrogramStream.fail())
105 {
106 m_electrogramStream >> gloCoord[0] >> gloCoord[1] >> gloCoord[2];
107 if (!m_electrogramStream.fail())
108 {
111 dim, i, gloCoord[0], gloCoord[1], gloCoord[2]);
112
113 m_electrogramPoints.push_back(vert);
114 ++i;
115 }
116 }
117
118 if (vComm->GetRank() == 0)
119 {
120 // Open output stream
121 m_outputStream.open(m_outputFile.c_str());
122 m_outputStream << "# Electrogram data for variables (:";
123
124 for (size_t i = 0; i < pFields.size(); ++i)
125 {
126 m_outputStream << m_session->GetVariable(i) << ",";
127 }
128
129 m_outputStream << ") at points:" << endl;
130
131 for (size_t i = 0; i < m_electrogramPoints.size(); ++i)
132 {
133 m_electrogramPoints[i]->GetCoords(gloCoord[0], gloCoord[1],
134 gloCoord[2]);
135
136 m_outputStream << "# \t" << i;
137 m_outputStream.width(8);
138 m_outputStream << gloCoord[0];
139 m_outputStream.width(8);
140 m_outputStream << gloCoord[1];
141 m_outputStream.width(8);
142 m_outputStream << gloCoord[2];
143 m_outputStream << endl;
144 }
145 }
146
147 // Compute the distance function for each electrogram point
148 const size_t nq = pFields[0]->GetNpoints();
149 const size_t npts = m_electrogramPoints.size();
150 NekDouble px, py, pz;
154
155 Array<OneD, NekDouble> oneOverR(nq, 0.0);
156
157 for (size_t i = 0; i < npts; ++i)
158 {
162
163 Array<OneD, NekDouble> x(nq, 0.0);
164 Array<OneD, NekDouble> y(nq, 0.0);
165 Array<OneD, NekDouble> z(nq, 0.0);
166
167 // Compute 1/R
168 m_electrogramPoints[i]->GetCoords(px, py, pz);
169 pFields[0]->GetCoords(x, y, z);
170
171 Vmath::Sadd(nq, -px, x, 1, x, 1);
172 Vmath::Sadd(nq, -py, y, 1, y, 1);
173 Vmath::Sadd(nq, -pz, z, 1, z, 1);
174 Vmath::Vvtvvtp(nq, x, 1, x, 1, y, 1, y, 1, oneOverR, 1);
175 Vmath::Vvtvp(nq, z, 1, z, 1, oneOverR, 1, oneOverR, 1);
176 Vmath::Vsqrt(nq, oneOverR, 1, oneOverR, 1);
177 Vmath::Sdiv(nq, 1.0, oneOverR, 1, oneOverR, 1);
178
179 // Compute grad 1/R
180 pFields[0]->PhysDeriv(oneOverR, m_grad_R_x[i], m_grad_R_y[i],
181 m_grad_R_z[i]);
182 }
183
184 // Compute electrogram point for initial condition
185 v_Update(pFields, time);
186}
187
188/**
189 *
190 */
193 const NekDouble &time)
194{
195 // Only output every m_outputFrequency.
196 if ((m_index++) % m_outputFrequency)
197 {
198 return;
199 }
200
201 const size_t nq = pFields[0]->GetNpoints();
202 const size_t npoints = m_electrogramPoints.size();
203 LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
204
205 size_t i = 0;
206 Array<OneD, NekDouble> e(npoints);
207
208 // Compute grad V
209 Array<OneD, NekDouble> grad_V_x(nq), grad_V_y(nq), grad_V_z(nq);
210 pFields[0]->PhysDeriv(pFields[0]->GetPhys(), grad_V_x, grad_V_y, grad_V_z);
211
212 for (i = 0; i < npoints; ++i)
213 {
214 // Multiply together
215 Array<OneD, NekDouble> output(nq);
216 Vmath::Vvtvvtp(nq, m_grad_R_x[i], 1, grad_V_x, 1, m_grad_R_y[i], 1,
217 grad_V_y, 1, output, 1);
218 Vmath::Vvtvp(nq, m_grad_R_z[i], 1, grad_V_z, 1, output, 1, output, 1);
219
220 e[i] = pFields[0]->Integral(output);
221 }
222
223 // Exchange history data
224 // This could be improved to reduce communication but works for now
225 vComm->AllReduce(e, LibUtilities::ReduceSum);
226
227 // Only the root process writes out electrogram data
228 if (vComm->GetRank() == 0)
229 {
230 m_outputStream.width(8);
231 m_outputStream << setprecision(6) << time;
232
233 // Write data values point by point
234 for (i = 0; i < m_electrogramPoints.size(); ++i)
235 {
236 m_outputStream.width(25);
237 m_outputStream << setprecision(16) << e[i];
238 }
239 m_outputStream << endl;
240 }
241}
242
243/**
244 *
245 */
248 [[maybe_unused]] const NekDouble &time)
249{
250 if (pFields[0]->GetComm()->GetRank() == 0)
251 {
252 m_outputStream.close();
253 }
254}
255
256/**
257 *
258 */
260{
261 return true;
262}
263} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
Array< OneD, Array< OneD, NekDouble > > m_grad_R_z
Gradient of the radius from each electrogram point in z-direction.
Array< OneD, Array< OneD, NekDouble > > m_grad_R_x
Gradient of the radius from each electrogram point in x-direction.
std::ofstream m_outputStream
Output file stream for electrogram data.
static SolverUtils::FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)
Creates an instance of this class.
std::string m_outputFile
Filename to output electrogram data to.
Array< OneD, Array< OneD, NekDouble > > m_grad_R_y
Gradient of the radius from each electrogram point in y-direction.
unsigned int m_outputFrequency
Number of timesteps between outputs.
~FilterElectrogram() override
Electrogram filter destructor.
FilterElectrogram(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)
Electrogram filter constructor.
bool v_IsTimeDependent() override
Filter is time-dependent and should be called at each time-step.
unsigned int m_index
Counts number of calls to update (number of timesteps)
SpatialDomains::PointGeomVector m_electrogramPoints
List of electrogram points.
void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
Compute extracellular potential at egm points at current time.
std::stringstream m_electrogramStream
Point coordinate input string.
static std::string className
Name of the class.
void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
Initialises the electrogram filter and open output file.
void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
Finalise the electrogram filter and close output file.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
SOLVER_UTILS_EXPORT std::string SetupOutput(const std::string ext, const ParamMap &pParams)
Definition: Filter.h:139
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:93
std::map< std::string, std::string > ParamMap
Definition: Filter.h:66
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
FilterFactory & GetFilterFactory()
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::vector< double > z(NPUPPER)
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
Definition: Vmath.hpp:154
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.hpp:439
STL namespace.