Skip to content

Commit eb6ffb8

Browse files
gateshiangateshian
committed
Implemented coronary outflow boundary condition
1 parent c27418f commit eb6ffb8

3 files changed

Lines changed: 356 additions & 0 deletions

File tree

‎FEBioFluid/FEBioFluid.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,7 @@ SOFTWARE.*/
5858
#include "FEFluidResistanceLoad.h"
5959
#include "FEFluidRCRBC.h"
6060
#include "FEFluidRCRLoad.h"
61+
#include "FEFluidCOBC.h"
6162
#include "FETangentialDamping.h"
6263
#include "FETangentialFlowStabilization.h"
6364
#include "FEBackFlowStabilization.h"
@@ -186,6 +187,7 @@ REGISTER_FECORE_CLASS(FEPrescribedFluidPressure , "fluid pressure");
186187
REGISTER_FECORE_CLASS(FEFluidRCBC , "fluid RC");
187188
REGISTER_FECORE_CLASS(FEFluidRCRBC , "fluid RCR");
188189
REGISTER_FECORE_CLASS(FEFluidResistanceBC , "fluid resistance");
190+
REGISTER_FECORE_CLASS(FEFluidCOBC , "fluid coronary outflow");
189191

190192
//-----------------------------------------------------------------------------
191193
// initial conditions

‎FEBioFluid/FEFluidCOBC.cpp‎

Lines changed: 256 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,256 @@
1+
/*This file is part of the FEBio source code and is licensed under the MIT license
2+
listed below.
3+
4+
See Copyright-FEBio.txt for details.
5+
6+
Copyright (c) 2021 University of Utah, The Trustees of Columbia University in
7+
the City of New York, and others.
8+
9+
Permission is hereby granted, free of charge, to any person obtaining a copy
10+
of this software and associated documentation files (the "Software"), to deal
11+
in the Software without restriction, including without limitation the rights
12+
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13+
copies of the Software, and to permit persons to whom the Software is
14+
furnished to do so, subject to the following conditions:
15+
16+
The above copyright notice and this permission notice shall be included in all
17+
copies or substantial portions of the Software.
18+
19+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20+
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22+
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24+
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25+
SOFTWARE.*/
26+
27+
#include "stdafx.h"
28+
#include "FEFluidCOBC.h"
29+
#include "FEFluid.h"
30+
#include "FEBioFluid.h"
31+
#include <FECore/FEAnalysis.h>
32+
#include <FECore/log.h>
33+
#include <FECore/FEModel.h>
34+
35+
//=============================================================================
36+
BEGIN_FECORE_CLASS(FEFluidCOBC, FEPrescribedSurface)
37+
ADD_PARAMETER(m_Ra, "Ra")->setLongName("arterial resistance")->setUnits("F.t/L^5");
38+
ADD_PARAMETER(m_Ca, "Ca")->setLongName("arterial compliance")->setUnits("L^5/F");
39+
ADD_PARAMETER(m_Ram,"Ra-micro")->setLongName("microvascular arterial resistance")->setUnits("F.t/L^5");
40+
ADD_PARAMETER(m_Cmi,"Cmi")->setLongName("myocardial compliance")->setUnits("L^5/F");
41+
ADD_PARAMETER(m_Rv, "Rv")->setLongName("ventricular resistance")->setUnits("F.t/L^5");
42+
ADD_PARAMETER(m_Rvm,"Rv-micro")->setLongName("microvascular ventricular resistance")->setUnits("F.t/L^5");
43+
ADD_PARAMETER(m_p0, "initial_pressure")->setUnits(UNIT_PRESSURE);
44+
// ADD_PARAMETER(m_pd, "pressure_offset")->setUnits(UNIT_PRESSURE);
45+
46+
ADD_PROPERTY(m_Pmi,"Pmi",FEProperty::Optional)->SetLongName("mycoardial pressure");
47+
ADD_PROPERTY(m_PRA,"PRA",FEProperty::Optional)->SetLongName("right atrium pressure");
48+
END_FECORE_CLASS();
49+
50+
//-----------------------------------------------------------------------------
51+
//! constructor
52+
FEFluidCOBC::FEFluidCOBC(FEModel* pfem) : FEPrescribedSurface(pfem), m_dofW(pfem)
53+
{
54+
m_Ra = m_Ram = m_Rv = m_Rvm = 0.0;
55+
m_Ca = m_Cmi = 0.0;
56+
m_pfluid = nullptr;
57+
m_psurf = nullptr;
58+
m_p0 = 0;
59+
m_pd = 0.0;
60+
m_e = 0.0;
61+
62+
}
63+
64+
//-----------------------------------------------------------------------------
65+
//! initialize
66+
//! TODO: Generalize to include the initial conditions
67+
bool FEFluidCOBC::Init()
68+
{
69+
m_dofW.AddVariable(FEBioFluid::GetVariableName(FEBioFluid::RELATIVE_FLUID_VELOCITY));
70+
m_dofEF = GetDOFIndex(FEBioFluid::GetVariableName(FEBioFluid::FLUID_DILATATION), 0);
71+
SetDOFList(m_dofEF);
72+
73+
if (FEPrescribedSurface::Init() == false) return false;
74+
75+
m_psurf = GetSurface();
76+
77+
// get fluid from first surface element
78+
// assuming the entire surface bounds the same fluid
79+
FESurfaceElement& el = m_psurf->Element(0);
80+
FEElement* pe = el.m_elem[0].pe;
81+
if (pe == nullptr) return false;
82+
83+
// get the material
84+
FEMaterial* pm = GetFEModel()->GetMaterial(pe->GetMatID());
85+
m_pfluid = pm->ExtractProperty<FEFluidMaterial>();
86+
if (m_pfluid == nullptr) return false;
87+
88+
m_pn = m_pp = m_ppp = m_p0;
89+
// m_pdn = m_pdp = m_pd;
90+
m_qn = m_qp = m_qpp = 0;
91+
m_tp = m_tpp = 0;
92+
93+
if (m_Pmi) { if (!m_Pmi->Init()) return false; }
94+
if (m_PRA) { if (!m_PRA->Init()) return false; }
95+
96+
return true;
97+
}
98+
99+
//-----------------------------------------------------------------------------
100+
//! Evaluate and prescribe the resistance pressure
101+
void FEFluidCOBC::UpdateDilatation()
102+
{
103+
// Check if we started a new time, if so, update variables
104+
FETimeInfo& timeInfo = GetFEModel()->GetTime();
105+
double time = timeInfo.currentTime;
106+
int iter = timeInfo.currentIteration;
107+
double dt = timeInfo.timeIncrement;
108+
if ((time > m_tp) && (iter == 0)) {
109+
m_ppp = m_pp;
110+
m_pp = m_pn;
111+
m_qpp = m_qp;
112+
m_qp = m_qn;
113+
m_pdp = m_pdn;
114+
m_tpp = m_tp;
115+
m_tp = time;
116+
}
117+
118+
// evaluate the flow rate at the current time
119+
m_qn = FlowRate();
120+
m_pdn = m_pd;
121+
122+
double Rve = m_Rv + m_Rvm;
123+
double Rae = m_Ram + Rve;
124+
double Rt = m_Ra + Rae;
125+
double tauvi = Rve*m_Cmi;
126+
double taua = m_Ram*m_Ca;
127+
double tauae = Rae*m_Ca;
128+
129+
double dtp = m_tp - m_tpp;
130+
131+
double c0 = tauvi*taua/(dt*dt);
132+
double c1 = m_Ra*c0;
133+
double c2 = (m_Ra*(tauae + tauvi)+m_Ram*tauvi)/dt;
134+
double a = c1 + c2 + Rt;
135+
double b = (c1+c2)*m_qp + m_Ra*tauvi*taua/dt/dtp*(m_qp-m_qpp);
136+
double c = c0 + (tauvi+tauae)/dt + 1;
137+
double d = c0*m_pp + tauvi*taua/dt/dtp*(m_pp - m_ppp) + (tauvi+tauae)/dt*m_pp;
138+
if (m_PRA) d += m_PRA->value(time);
139+
if (m_Pmi) d += tauvi*m_Pmi->derive(time);
140+
141+
// calculate the RCR pressure
142+
m_pn = (a*m_qn + d - b)/c;
143+
144+
// calculate the dilatation
145+
m_e = 0.0;
146+
bool good = m_pfluid->Dilatation(0, m_pn, m_e);
147+
assert(good);
148+
}
149+
150+
void FEFluidCOBC::UpdateModel() { Update(); }
151+
152+
void FEFluidCOBC::Update()
153+
{
154+
UpdateDilatation();
155+
156+
// the base class handles mapping the values to the nodal dofs
157+
FEPrescribedSurface::Update();
158+
159+
// TODO: Is this necessary?
160+
GetFEModel()->SetMeshUpdateFlag(true);
161+
}
162+
163+
//-----------------------------------------------------------------------------
164+
//! evaluate the flow rate across this surface at current time
165+
double FEFluidCOBC::FlowRate()
166+
{
167+
double Q = 0;
168+
169+
const FETimeInfo& tp = GetTimeInfo();
170+
171+
vec3d rt[FEElement::MAX_NODES];
172+
vec3d vt[FEElement::MAX_NODES];
173+
174+
for (int iel=0; iel<m_psurf->Elements(); ++iel)
175+
{
176+
FESurfaceElement& el = m_psurf->Element(iel);
177+
178+
// nr integration points
179+
int nint = el.GaussPoints();
180+
181+
// nr of element nodes
182+
int neln = el.Nodes();
183+
184+
// nodal coordinates
185+
for (int i=0; i<neln; ++i) {
186+
FENode& node = m_psurf->Node(el.m_lnode[i]);
187+
rt[i] = node.m_rt;
188+
vt[i] = node.get_vec3d(m_dofW[0], m_dofW[1], m_dofW[2]);
189+
}
190+
191+
double* Nr, *Ns;
192+
double* N;
193+
double* w = el.GaussWeights();
194+
195+
vec3d dxr, dxs, v;
196+
197+
// repeat over integration points
198+
for (int n=0; n<nint; ++n)
199+
{
200+
N = el.H(n);
201+
Nr = el.Gr(n);
202+
Ns = el.Gs(n);
203+
204+
// calculate the velocity and tangent vectors at integration point
205+
dxr = dxs = v = vec3d(0,0,0);
206+
for (int i=0; i<neln; ++i)
207+
{
208+
v += vt[i]*N[i];
209+
dxr += rt[i]*Nr[i];
210+
dxs += rt[i]*Ns[i];
211+
}
212+
213+
vec3d normal = dxr ^ dxs;
214+
double q = normal*v;
215+
Q += q*w[n];
216+
}
217+
}
218+
219+
return Q;
220+
}
221+
222+
//-----------------------------------------------------------------------------
223+
void FEFluidCOBC::PrepStep(std::vector<double>& ui, bool brel)
224+
{
225+
UpdateDilatation();
226+
FEPrescribedSurface::PrepStep(ui, brel);
227+
}
228+
229+
//-----------------------------------------------------------------------------
230+
void FEFluidCOBC::GetNodalValues(int nodelid, std::vector<double>& val)
231+
{
232+
val[0] = m_e;
233+
FENode& node = GetMesh().Node(m_nodeList[nodelid]);
234+
node.set(m_dofEF, m_e);
235+
}
236+
237+
//-----------------------------------------------------------------------------
238+
// copy data from another class
239+
void FEFluidCOBC::CopyFrom(FEBoundaryCondition* pbc)
240+
{
241+
// TODO: implement this
242+
assert(false);
243+
}
244+
245+
//-----------------------------------------------------------------------------
246+
//! serialization
247+
void FEFluidCOBC::Serialize(DumpStream& ar)
248+
{
249+
FEPrescribedSurface::Serialize(ar);
250+
ar & m_pn & m_pp & m_pp & m_qn & m_qp & m_qpp & m_tp & m_tpp & m_e;
251+
if (ar.IsShallow()) return;
252+
ar & m_pfluid;
253+
ar & m_dofW & m_dofEF;
254+
ar & m_psurf;
255+
ar & m_Pmi & m_PRA;
256+
}

‎FEBioFluid/FEFluidCOBC.h‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
/*This file is part of the FEBio source code and is licensed under the MIT license
2+
listed below.
3+
4+
See Copyright-FEBio.txt for details.
5+
6+
Copyright (c) 2021 University of Utah, The Trustees of Columbia University in
7+
the City of New York, and others.
8+
9+
Permission is hereby granted, free of charge, to any person obtaining a copy
10+
of this software and associated documentation files (the "Software"), to deal
11+
in the Software without restriction, including without limitation the rights
12+
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13+
copies of the Software, and to permit persons to whom the Software is
14+
furnished to do so, subject to the following conditions:
15+
16+
The above copyright notice and this permission notice shall be included in all
17+
copies or substantial portions of the Software.
18+
19+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20+
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22+
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24+
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25+
SOFTWARE.*/
26+
27+
#pragma once
28+
#include <FECore/FEPrescribedBC.h>
29+
#include <FECore/FEFunction1D.h>
30+
#include "FEFluidMaterial.h"
31+
32+
//-----------------------------------------------------------------------------
33+
//! FEFluidCOBC is a fluid surface load that implements a coronary outflow boundary condition
34+
//!
35+
class FEBIOFLUID_API FEFluidCOBC : public FEPrescribedSurface
36+
{
37+
public:
38+
//! constructor
39+
FEFluidCOBC(FEModel* pfem);
40+
41+
//! set the dilatation
42+
void Update() override;
43+
void UpdateModel() override;
44+
45+
//! evaluate flow rate
46+
double FlowRate();
47+
48+
//! initialize
49+
bool Init() override;
50+
51+
//! serialization
52+
void Serialize(DumpStream& ar) override;
53+
54+
public:
55+
void PrepStep(std::vector<double>& ui, bool brel) override;
56+
57+
// return the value for node i, dof j
58+
void GetNodalValues(int nodelid, std::vector<double>& val) override;
59+
60+
// copy data from another class
61+
void CopyFrom(FEBoundaryCondition* pbc) override;
62+
63+
private:
64+
//! set the dilatation
65+
void UpdateDilatation();
66+
67+
private:
68+
double m_Ra; //!< arterial resistance
69+
double m_Ca; //!< arterial compliance
70+
double m_Ram; //!< microvascular arterial resistance
71+
double m_Cmi; //!< myocardial compliance
72+
double m_Rv; //!< ventricular resistance
73+
double m_Rvm; //!< microvascular ventricular resistance
74+
double m_p0; //!< initial fluid pressure
75+
FEFunction1D* m_Pmi; //!< myocardial pressure
76+
FEFunction1D* m_PRA; //!< right atrium pressure
77+
double m_pd; //!< downstream pressure
78+
79+
private:
80+
double m_pn; //!< fluid pressure at current time point
81+
double m_pp; //!< fluid pressure at previous time point
82+
double m_ppp; //!< fluid pressure at 2nd-previous time point
83+
double m_qn; //!< flow rate at current time point
84+
double m_qp; //!< flow rate at previous time point
85+
double m_qpp; //!< flow rate at 2nd-previous time point
86+
double m_pdn; //!< downstream fluid pressure at current time point
87+
double m_pdp; //!< downstream fluid pressure at previous time point
88+
double m_tp; //!< previous time
89+
double m_tpp; //!< 2nd-previous time
90+
double m_e;
91+
FEFluidMaterial* m_pfluid; //!< pointer to fluid
92+
FESurface* m_psurf; //!< pointer to surface
93+
94+
FEDofList m_dofW;
95+
int m_dofEF;
96+
97+
DECLARE_FECORE_CLASS();
98+
};

0 commit comments

Comments
 (0)