44 "Beta law pressure area relationship for the arterial system");
60 gamma * dAUdx /
sqrt(
A);
127 const std::string &type)
138 if (type ==
"Bifurcation")
144 for (
int i = 0; i < 3; ++i)
146 GetC(c[i],
beta[i], Au[i], A0[i], alpha[i]);
149 k = c[0] * Au[1] * c[2] + Au[0] * c[2] * c[1] + Au[2] * c[0] * c[1];
150 K[0] = (c[0] - uu[0]) * k;
151 K[1] = (c[1] + uu[1]) * k;
152 K[2] = (c[2] + uu[2]) * k;
155 (-c[1] * uu[0] * c[2] * Au[0] +
156 Au[2] * c[1] * c[0] * c[0] +
157 Au[1] * c[0] * c[0] * c[2]) /
159 invJ.SetValue(0, 1, Au[1] * (c[1] - uu[1]) * c[0] * c[2] / K[0]);
160 invJ.SetValue(0, 2, Au[2] * (c[2] - uu[2]) * c[0] * c[1] / K[0]);
161 invJ.SetValue(0, 3, c[0] * c[1] * c[2] / K[0]);
162 invJ.SetValue(0, 4, -0.5 * c[0] * Au[1] * c[2] / K[0]);
163 invJ.SetValue(0, 5, -0.5 * Au[2] * c[0] * c[1] / K[0]);
165 invJ.SetValue(1, 0, Au[0] * (c[0] + uu[0]) * c[1] * c[2] / K[1]);
167 (c[0] * uu[1] * c[2] * Au[1] +
168 Au[2] * c[0] * c[1] * c[1] +
169 c[2] * c[1] * c[1] * Au[0]) /
171 invJ.SetValue(1, 2, -Au[2] * (c[2] - uu[2]) * c[0] * c[1] / K[1]);
172 invJ.SetValue(1, 3, -c[0] * c[1] * c[2] / K[1]);
173 invJ.SetValue(1, 4, -0.5 * (c[0] * Au[2] + Au[0] * c[2]) * c[1] / K[1]);
174 invJ.SetValue(1, 5, 0.5 * Au[2] * c[0] * c[1] / K[1]);
176 invJ.SetValue(2, 0, Au[0] * (c[0] + uu[0]) * c[1] * c[2] / K[2]);
177 invJ.SetValue(2, 1, -Au[1] * (c[1] - uu[1]) * c[0] * c[2] / K[2]);
179 (c[0] * c[1] * uu[2] * Au[2] +
180 c[0] * Au[1] * c[2] * c[2] +
181 c[1] * c[2] * c[2] * Au[0]) /
183 invJ.SetValue(2, 3, -c[0] * c[1] * c[2] / K[2]);
184 invJ.SetValue(2, 4, 0.5 * c[0] * Au[1] * c[2] / K[2]);
185 invJ.SetValue(2, 5, -0.5 * (Au[1] * c[0] + c[1] * Au[0]) * c[2] / K[2]);
189 (Au[0] * c[2] * c[1] - uu[0] * c[2] * Au[1] -
190 uu[0] * c[1] * Au[2]) /
192 invJ.SetValue(3, 1, -Au[0] * Au[1] * (c[1] - uu[1]) * c[2] / K[0]);
193 invJ.SetValue(3, 2, -Au[0] * Au[2] * (c[2] - uu[2]) * c[1] / K[0]);
194 invJ.SetValue(3, 3, -Au[0] * c[2] * c[1] / K[0]);
195 invJ.SetValue(3, 4, 0.5 * Au[0] * Au[1] * c[2] / K[0]);
196 invJ.SetValue(3, 5, 0.5 * Au[0] * c[1] * Au[2] / K[0]);
198 invJ.SetValue(4, 0, Au[0] * Au[1] * (c[0] + uu[0]) * c[2] / K[1]);
201 (c[0] * Au[1] * c[2] + c[0] * uu[1] * Au[2] +
202 c[2] * uu[1] * Au[0]) /
204 invJ.SetValue(4, 2, -Au[2] * Au[1] * (c[2] - uu[2]) * c[0] / K[1]);
205 invJ.SetValue(4, 3, -c[0] * Au[1] * c[2] / K[1]);
207 -0.5 * Au[1] * (c[0] * Au[2] + Au[0] * c[2]) / K[1]);
208 invJ.SetValue(4, 5, 0.5 * Au[2] * Au[1] * c[0] / K[1]);
210 invJ.SetValue(5, 0, Au[0] * Au[2] * (c[0] + uu[0]) * c[1] / K[2]);
211 invJ.SetValue(5, 1, -Au[2] * Au[1] * (c[1] - uu[1]) * c[0] / K[2]);
214 (Au[2] * c[0] * c[1] + c[0] * uu[2] * Au[1] +
215 c[1] * uu[2] * Au[0]) /
217 invJ.SetValue(5, 3, -Au[2] * c[0] * c[1] / K[2]);
218 invJ.SetValue(5, 4, 0.5 * Au[2] * Au[1] * c[0] / K[2]);
220 -0.5 * Au[2] * (Au[1] * c[0] + c[1] * Au[0]) / K[2]);
222 else if (type ==
"Merge")
228 for (
int i = 0; i < 3; ++i)
230 GetC(c[i],
beta[i], Au[i], A0[i], alpha[i]);
233 k = c[0] * Au[1] * c[2] + Au[0] * c[2] * c[1] + Au[2] * c[0] * c[1];
234 K[0] = (c[0] - uu[0]) * k;
235 K[1] = (c[1] + uu[1]) * k;
236 K[2] = (c[2] + uu[2]) * k;
239 (-c[1] * uu[0] * c[2] * Au[0] +
240 Au[2] * c[1] * c[0] * c[0] +
241 Au[1] * c[0] * c[0] * c[2]) /
243 invJ.SetValue(0, 1, Au[1] * (c[1] - uu[1]) * c[0] * c[2] / K[0]);
244 invJ.SetValue(0, 2, Au[2] * (c[2] - uu[2]) * c[0] * c[1] / K[0]);
245 invJ.SetValue(0, 3, c[0] * c[1] * c[2] / K[0]);
246 invJ.SetValue(0, 4, -0.5 * c[0] * Au[1] * c[2] / K[0]);
247 invJ.SetValue(0, 5, -0.5 * Au[2] * c[0] * c[1] / K[0]);
249 invJ.SetValue(1, 0, Au[0] * (c[0] + uu[0]) * c[1] * c[2] / K[1]);
251 (c[0] * uu[1] * c[2] * Au[1] +
252 Au[2] * c[0] * c[1] * c[1] +
253 c[2] * c[1] * c[1] * Au[0]) /
255 invJ.SetValue(1, 2, -Au[2] * (c[2] - uu[2]) * c[0] * c[1] / K[1]);
256 invJ.SetValue(1, 3, -c[0] * c[1] * c[2] / K[1]);
257 invJ.SetValue(1, 4, -0.5 * (c[0] * Au[2] + Au[0] * c[2]) * c[1] / K[1]);
258 invJ.SetValue(1, 5, 0.5 * Au[2] * c[0] * c[1] / K[1]);
260 invJ.SetValue(2, 0, Au[0] * (c[0] - uu[0]) * c[1] * c[2] / K[2]);
261 invJ.SetValue(2, 1, -Au[1] * (c[1] + uu[1]) * c[0] * c[2] / K[2]);
263 -(c[0] * uu[2] * c[1] * Au[2] -
264 Au[1] * c[0] * c[2] * c[2] -
265 c[1] * c[2] * c[2] * Au[0]) /
267 invJ.SetValue(2, 3, -c[0] * c[1] * c[2] / K[2]);
268 invJ.SetValue(2, 4, -0.5 * Au[1] * c[0] * c[2] / K[2]);
269 invJ.SetValue(2, 5, 0.5 * (Au[1] * c[0] + Au[0] * c[1]) * c[2] / K[2]);
273 (Au[0] * c[2] * c[1] + uu[0] * c[2] * Au[1] +
274 uu[0] * c[1] * Au[2]) /
276 invJ.SetValue(3, 1, Au[0] * Au[1] * (c[1] + uu[1]) * c[2] / K[0]);
278 invJ.SetValue(3, 2, -Au[0] * Au[2] * (c[2] - uu[2]) * c[1] / K[0]);
279 invJ.SetValue(3, 3, -Au[0] * c[2] * c[1] / K[0]);
280 invJ.SetValue(3, 4, 0.5 * Au[0] * Au[1] * c[2] / K[0]);
281 invJ.SetValue(3, 5, 0.5 * Au[0] * c[1] * Au[2] / K[0]);
283 invJ.SetValue(4, 0, Au[0] * Au[1] * (c[0] + uu[0]) * c[2] / K[1]);
286 (c[0] * Au[1] * c[2] + c[0] * uu[1] * Au[2] +
287 c[2] * uu[1] * Au[0]) /
289 invJ.SetValue(4, 2, -Au[2] * Au[1] * (c[2] - uu[2]) * c[0] / K[1]);
290 invJ.SetValue(4, 3, -c[0] * Au[1] * c[2] / K[1]);
292 -0.5 * Au[1] * (c[0] * Au[2] + Au[0] * c[2]) / K[1]);
293 invJ.SetValue(4, 5, 0.5 * Au[2] * Au[1] * c[0] / K[1]);
295 invJ.SetValue(5, 0, Au[0] * Au[2] * (c[0] + uu[0]) * c[1] / K[2]);
296 invJ.SetValue(5, 1, -Au[2] * Au[1] * (c[1] - uu[1]) * c[0] / K[2]);
299 (Au[2] * c[0] * c[1] + c[0] * uu[2] * Au[1] +
300 c[1] * uu[2] * Au[0]) /
302 invJ.SetValue(5, 3, -Au[2] * c[0] * c[1] / K[2]);
303 invJ.SetValue(5, 4, 0.5 * Au[2] * Au[1] * c[0] / K[2]);
305 -0.5 * Au[2] * (Au[1] * c[0] + c[1] * Au[0]) / K[2]);
307 else if (type ==
"Interface")
313 for (
int i = 0; i < 2; ++i)
315 GetC(c[i],
beta[i], Au[i], A0[i], alpha[i]);
318 k = (c[0] * Au[1] + Au[0] * c[1]);
319 K[0] = (c[0] + uu[0]) * k;
320 K[1] = (c[1] + uu[1]) * k;
323 (Au[1] * c[0] * c[0] - c[1] * uu[0] * Au[0]) / K[0]);
324 invJ.SetValue(0, 1, Au[1] * (c[1] - uu[1]) * c[0] / K[0]);
325 invJ.SetValue(0, 2, c[0] * c[1] / K[0]);
326 invJ.SetValue(0, 3, -0.5 * c[0] * Au[1] / K[0]);
328 invJ.SetValue(1, 0, Au[0] * (c[0] + uu[0]) * c[1] / K[1]);
330 (c[0] * uu[1] * Au[1] + c[1] * c[1] * Au[0]) / K[1]);
331 invJ.SetValue(1, 2, -c[0] * c[1] / K[1]);
332 invJ.SetValue(1, 3, -0.5 * Au[0] * c[1] / K[1]);
334 invJ.SetValue(2, 0, Au[0] * (Au[0] * c[1] - uu[0] * Au[1]) / K[0]);
335 invJ.SetValue(2, 1, -Au[0] * Au[1] * (c[1] - uu[1]) / K[0]);
336 invJ.SetValue(2, 2, -Au[0] * c[1] / K[0]);
337 invJ.SetValue(2, 3, 0.5 * Au[1] * Au[0] / K[0]);
339 invJ.SetValue(3, 0, Au[0] * Au[1] * (c[0] + uu[0]) / K[1]);
340 invJ.SetValue(3, 1, -Au[1] * (c[0] * Au[1] + uu[1] * Au[0]) / K[1]);
341 invJ.SetValue(3, 2, -c[0] * Au[1] / K[1]);
342 invJ.SetValue(3, 3, -0.5 * Au[1] * Au[0] / K[1]);
BetaPressureArea(Array< OneD, MultiRegions::ExpListSharedPtr > pVessel, const LibUtilities::SessionReaderSharedPtr pSession)
void v_GetPressure(NekDouble &P, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &dAUdx, const NekDouble &gamma=0, const NekDouble &alpha=0.5) override
void v_GetCharIntegral(NekDouble &I, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5) override
static PulseWavePressureAreaSharedPtr create(Array< OneD, MultiRegions::ExpListSharedPtr > &pVessel, const LibUtilities::SessionReaderSharedPtr &pSession)
static std::string className
void v_GetW2(NekDouble &W2, const NekDouble &u, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5) override
void v_GetUFromChars(NekDouble &u, const NekDouble &W1, const NekDouble &W2) override
void v_GetC(NekDouble &c, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5) override
void v_GetJacobianInverse(NekMatrix< NekDouble > &invJ, const Array< OneD, NekDouble > &Au, const Array< OneD, NekDouble > &uu, const Array< OneD, NekDouble > &beta, const Array< OneD, NekDouble > &A0, const Array< OneD, NekDouble > &alpha, const std::string &type) override
void v_GetW1(NekDouble &W1, const NekDouble &u, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5) override
void v_GetAFromChars(NekDouble &A, const NekDouble &W1, const NekDouble &W2, const NekDouble &beta, const NekDouble &A0, const NekDouble &alpha=0.5) override
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
void GetCharIntegral(NekDouble &I, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
void GetC(NekDouble &c, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
@ P
Monomial polynomials .
PressureAreaFactory & GetPressureAreaFactory()
scalarT< T > sqrt(scalarT< T > in)