Skip to content

Commit

Permalink
Option to weigh active stress field in FrSI
Browse files Browse the repository at this point in the history
  • Loading branch information
Marc Hirschvogel committed Sep 25, 2024
1 parent 474f354 commit f181344
Show file tree
Hide file tree
Showing 5 changed files with 42 additions and 15 deletions.
12 changes: 9 additions & 3 deletions src/ambit_fe/boundaryconditions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -485,7 +485,7 @@ def robin_bcs(self, bcdict, u, v, ds_, u_pre=None, wel=None, F=None):


# set membrane surface BCs
def membranesurf_bcs(self, bcdict, u, v, a, ds_, ivar=None, wallfields=[]):
def membranesurf_bcs(self, bcdict, u, v, a, ds_, ivar=None, wallfields=[], actweights=[]):

w, idmem, bstress, bstrainenergy, bintpower = ufl.as_ufl(0), [], [], [], []

Expand Down Expand Up @@ -515,12 +515,18 @@ def membranesurf_bcs(self, bcdict, u, v, a, ds_, ivar=None, wallfields=[]):
else:
wallfield = None

# field for active stress weighting
if bool(actweights):
actweight = actweights[mi]
else:
actweight = None

for i in range(len(m['id'])):

idmem.append(m['id'][i])

w += self.vf.deltaW_ext_membrane(self.ki.F(u), self.ki.Fdot(v), a, m['params'], ds_[dind](m['id'][i]), ivar=ivar, fibfnc=self.ff, wallfield=wallfield, fcts=fcts)
bstr, bse, bip = self.vf.deltaW_ext_membrane(self.ki.F(u), self.ki.Fdot(v), a, m['params'], ds_[dind](m['id'][i]), ivar=ivar, fibfnc=self.ff, wallfield=wallfield, fcts=fcts, returnquantity='stress_energy_power')
w += self.vf.deltaW_ext_membrane(self.ki.F(u), self.ki.Fdot(v), a, m['params'], ds_[dind](m['id'][i]), ivar=ivar, fibfnc=self.ff, wallfield=wallfield, actweight=actweight, fcts=fcts)
bstr, bse, bip = self.vf.deltaW_ext_membrane(self.ki.F(u), self.ki.Fdot(v), a, m['params'], ds_[dind](m['id'][i]), ivar=ivar, fibfnc=self.ff, wallfield=wallfield, actweight=actweight, fcts=fcts, returnquantity='stress_energy_power')
bstress.append(bstr)
bstrainenergy.append(bse)
bintpower.append(bip)
Expand Down
18 changes: 13 additions & 5 deletions src/ambit_fe/coupling/fluid_constraint_main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -193,14 +193,22 @@ def set_variational_forms(self):
ivar_old_ = {'tau_a' : self.coupfuncs_old[-1]}
ivar_mid_ = {'tau_a' : self.coupfuncs_mid[-1]}

if 'weight' in self.coupling_params['multiplier_physics'][n].keys():
# active stress weighting for reduced solid
wact_func = fem.Function(self.pbf.V_scalar)
self.pbf.io.readfunction(wact_func, self.coupling_params['multiplier_physics'][n]['weight'])
self.pbf.actweights.append(wact_func)
else:
self.pbf.actweights.append(None)

# add internal active stress power to fluid rhs contributions
self.power_coupling += self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.ufluid), self.pbf.ki.Fdot(self.pbf.v), None, params_, ds_p, ivar=ivar_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], returnquantity='active_stress_power')
self.power_coupling_old += self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.uf_old), self.pbf.ki.Fdot(self.pbf.v_old), None, params_, ds_p, ivar=ivar_old_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], returnquantity='active_stress_power')
self.power_coupling_mid += self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.ufluid_mid), self.pbf.ki.Fdot(self.pbf.vel_mid), None, params_, ds_p, ivar=ivar_mid_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], returnquantity='active_stress_power')
self.power_coupling += self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.ufluid), self.pbf.ki.Fdot(self.pbf.v), None, params_, ds_p, ivar=ivar_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], actweight=self.pbf.actweights[-1], returnquantity='active_stress_power')
self.power_coupling_old += self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.uf_old), self.pbf.ki.Fdot(self.pbf.v_old), None, params_, ds_p, ivar=ivar_old_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], actweight=self.pbf.actweights[-1], returnquantity='active_stress_power')
self.power_coupling_mid += self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.ufluid_mid), self.pbf.ki.Fdot(self.pbf.vel_mid), None, params_, ds_p, ivar=ivar_mid_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], actweight=self.pbf.actweights[-1], returnquantity='active_stress_power')

# derivative w.r.t. multiplier
df_ += self.pbf.timefac * self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.ufluid), self.pbf.ki.Fdot(self.pbf.v), None, params_, ds_p, ivar=ivar_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], returnquantity='active_stress_power_deriv')
df_mid_ += self.pbf.timefac * self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.ufluid_mid), self.pbf.ki.Fdot(self.pbf.vel_mid), None, params_, ds_p, ivar=ivar_mid_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], returnquantity='active_stress_power_deriv')
df_ += self.pbf.timefac * self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.ufluid), self.pbf.ki.Fdot(self.pbf.v), None, params_, ds_p, ivar=ivar_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], actweight=self.pbf.actweights[-1], returnquantity='active_stress_power_deriv')
df_mid_ += self.pbf.timefac * self.pbf.vf.deltaW_ext_membrane(self.pbf.ki.F(self.pbf.ufluid_mid), self.pbf.ki.Fdot(self.pbf.vel_mid), None, params_, ds_p, ivar=ivar_mid_, fibfnc=self.pbf.fib_func, wallfield=self.pbf.wallfields[n], actweight=self.pbf.actweights[-1], returnquantity='active_stress_power_deriv')

elif self.coupling_params['multiplier_physics'][n]['type'] == 'valve_viscosity':

Expand Down
16 changes: 12 additions & 4 deletions src/ambit_fe/fluid/fluid_main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -588,7 +588,7 @@ def set_variational_forms(self):

self.internalvars['tau_a'], self.internalvars_old['tau_a'], self.internalvars_mid['tau_a'] = self.tau_a, self.tau_a_old, self.timefac*self.tau_a + (1.-self.timefac)*self.tau_a_old

self.actstress, self.act_curve, self.wallfields = [], [], []
self.actstress, self.act_curve, self.wallfields, self.actweights = [], [], [], []
for nm in range(len(self.bc_dict['membrane'])):

if 'active_stress' in self.bc_dict['membrane'][nm]['params'].keys():
Expand All @@ -610,6 +610,14 @@ def set_variational_forms(self):
else:
raise NameError("Unknown active stress type for membrane!")

if 'weight' in self.bc_dict['membrane'][nm]['params']['active_stress'].keys():
# active stress weighting for reduced solid
wact_func = fem.Function(self.V_scalar)
self.io.readfunction(wact_func, self.bc_dict['membrane'][nm]['params']['active_stress']['weight'])
self.actweights.append(wact_func)
else:
self.actweights.append(None)

if 'field' in self.bc_dict['membrane'][nm]['params']['h0'].keys():
# wall thickness field for reduced solid
h0_func = fem.Function(self.V_scalar)
Expand All @@ -618,9 +626,9 @@ def set_variational_forms(self):
else:
self.wallfields.append(None)

w_membrane, self.idmem, self.bstress, self.bstrainenergy, self.bintpower = self.bc.membranesurf_bcs(self.bc_dict['membrane'], self.ufluid, self.v, self.acc, self.bmeasures, ivar=self.internalvars, wallfields=self.wallfields)
w_membrane_old, _, _, _, _ = self.bc.membranesurf_bcs(self.bc_dict['membrane'], self.uf_old, self.v_old, self.a_old, self.bmeasures, ivar=self.internalvars_old, wallfields=self.wallfields)
w_membrane_mid, _, _, _, _ = self.bc.membranesurf_bcs(self.bc_dict['membrane'], self.ufluid_mid, self.vel_mid, self.acc_mid, self.bmeasures, ivar=self.internalvars_mid, wallfields=self.wallfields)
w_membrane, self.idmem, self.bstress, self.bstrainenergy, self.bintpower = self.bc.membranesurf_bcs(self.bc_dict['membrane'], self.ufluid, self.v, self.acc, self.bmeasures, ivar=self.internalvars, wallfields=self.wallfields, actweights=self.actweights)
w_membrane_old, _, _, _, _ = self.bc.membranesurf_bcs(self.bc_dict['membrane'], self.uf_old, self.v_old, self.a_old, self.bmeasures, ivar=self.internalvars_old, wallfields=self.wallfields, actweights=self.actweights)
w_membrane_mid, _, _, _, _ = self.bc.membranesurf_bcs(self.bc_dict['membrane'], self.ufluid_mid, self.vel_mid, self.acc_mid, self.bmeasures, ivar=self.internalvars_mid, wallfields=self.wallfields, actweights=self.actweights)

w_neumann_prestr, self.deltaW_prestr_kin, self.deltaW_prestr_int, self.deltaW_p_prestr = ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0), []
if self.prestress_initial or self.prestress_initial_only:
Expand Down
8 changes: 6 additions & 2 deletions src/ambit_fe/variationalform.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@ def deltaW_ext_robin_dashpot_normal_cross(self, v, c_c, dboundary):
# TeX: h_0\int\limits_{\Gamma_0} \boldsymbol{P}(\boldsymbol{u},\boldsymbol{v}(\boldsymbol{u})) : \boldsymbol{\nabla}_{\tilde{\boldsymbol{X}}}\delta\boldsymbol{u}\,\mathrm{d}A
# for fluid mechanics, contribution to virtual power is:
# TeX: h_0\int\limits_{\Gamma_0} \boldsymbol{P}(\boldsymbol{u}_{\mathrm{f}}(\boldsymbol{v}),\boldsymbol{v}) : \boldsymbol{\nabla}_{\tilde{\boldsymbol{X}}}\delta\boldsymbol{v}\,\mathrm{d}A
def deltaW_ext_membrane(self, F, Fdot, a, params, dboundary, ivar=None, fibfnc=None, wallfield=None, fcts=None, returnquantity='weakform'):
def deltaW_ext_membrane(self, F, Fdot, a, params, dboundary, ivar=None, fibfnc=None, wallfield=None, actweight=None, fcts=None, returnquantity='weakform'):

C = F.T*F

Expand All @@ -160,6 +160,10 @@ def deltaW_ext_membrane(self, F, Fdot, a, params, dboundary, ivar=None, fibfnc=N
dS_act = self.I
else:
ValueError("Unknown ative stress dir!")
if actweight is not None:
w_act = actweight
else:
w_act = 1.0

# wall thickness - can be constant or a field
wall_thickness = params['h0']
Expand Down Expand Up @@ -257,7 +261,7 @@ def deltaW_ext_membrane(self, F, Fdot, a, params, dboundary, ivar=None, fibfnc=N

# add active stress
if active is not None:
S += S_act
S += w_act * S_act

# 1st PK stress P = FS
P = Fmod * S
Expand Down
3 changes: 2 additions & 1 deletion tests/test_frsi_blocks_active.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

"""
FrSI test case of three blocks with an active membrane surface, paritioned POD space per block that allows in-plane contraction of the membrane
tests also active stress weight for first block (though using only the first partition field, hence multiply by 1)
"""

import ambit_fe
Expand Down Expand Up @@ -98,7 +99,7 @@ def tc2(self, t): # active stress from file - precomputed from the ODE, should y

BC_DICT_ALE = { 'dirichlet' : [{'id' : [5,11,17], 'dir' : 'all', 'val' : 0.}] } # bottom

BC_DICT_FLUID = { 'membrane' : [{'id' : [2], 'params' : {'model' : 'membrane', 'a_0' : 1.0, 'b_0' : 6.0, 'eta' : 0.01, 'rho0' : 1e-6, 'h0' : {'val' : 1.0}, 'active_stress' : {'type' : 'ode', 'dir' : 'cl', 'sigma0' : 10., 'alpha_max' : 10.0, 'alpha_min' : -30.0, 'activation_curve' : 1, 'omega' : 0.667, 'iota' : 0.333, 'gamma' : 0.0}}},
BC_DICT_FLUID = { 'membrane' : [{'id' : [2], 'params' : {'model' : 'membrane', 'a_0' : 1.0, 'b_0' : 6.0, 'eta' : 0.01, 'rho0' : 1e-6, 'h0' : {'val' : 1.0}, 'active_stress' : {'type' : 'ode', 'dir' : 'cl', 'sigma0' : 10., 'alpha_max' : 10.0, 'alpha_min' : -30.0, 'activation_curve' : 1, 'omega' : 0.667, 'iota' : 0.333, 'gamma' : 0.0, 'weight' : basepath+'/input/blocks3_part-1.txt'}}},
{'id' : [8], 'params' : {'model' : 'membrane', 'a_0' : 1.0, 'b_0' : 6.0, 'eta' : 0.01, 'rho0' : 1e-6, 'h0' : {'val' : 1.0}, 'active_stress' : {'type' : 'prescribed', 'dir' : 'cl', 'prescribed_curve' : 2, 'omega' : 0.667, 'iota' : 0.333, 'gamma' : 0.0}}},
{'id' : [14], 'params' : {'model' : 'membrane', 'a_0' : 1.0, 'b_0' : 6.0, 'eta' : 0.01, 'rho0' : 1e-6, 'h0' : {'val' : 1.0}, 'active_stress' : {'type' : 'ode', 'dir' : 'iso', 'sigma0' : 100., 'alpha_max' : 10.0, 'alpha_min' : -30.0, 'activation_curve' : 1}}}] }

Expand Down

0 comments on commit f181344

Please sign in to comment.