# This file is part of `dolfin_navier_scipy`.
#
# `dolfin_navier_scipy` is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# `dolfin_navier_scipy` is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with `dolfin_navier_scipy`.
# If not, see <http://www.gnu.org/licenses/>.
import logging
import dolfin
import dolfin_navier_scipy.dolfin_to_sparrays as dts
import numpy as np
import json
__all__ = ['get_sysmats',
'drivcav_fems',
'cyl_fems',
'gen_bccont_fems',
'gen_bccont_fems_3D',
'cyl3D_fems',
'get_bcinds',
]
[docs]def get_sysmats(problem='gen_bccont', scheme=None, ppin=None,
Re=None, nu=None, charvel=1., gradvsymmtrc=True,
bccontrol=False, mergerhs=False,
onlymesh=False, meshparams={}):
""" retrieve the system matrices for stokes flow
Parameters
----------
problem : {'drivencavity', 'cylinderwake', 'gen_bccont', 'cylinder_rot'}
problem class
N : int
mesh parameter
nu : real, optional
kinematic viscosity, is set to `L/Re` if `Re` is provided
Re : real, optional
Reynoldsnumber, is set to `L/nu` if `nu` is provided
bccontrol : boolean, optional
whether to consider boundary control via penalized Robin \
defaults to `False`
mergerhs : boolean, optional
whether to merge the actual rhs and the contribution from the \
boundary conditions into one rhs, defaults to `False`
onlymesh : boolean, optional
whether to only return `femp`, containing the mesh and FEM spaces, \
defaults to `False`
Returns
-------
femp : dict
with the keys:
* `V`: FEM space of the velocity
* `Q`: FEM space of the pressure
* `diribcs`: list of the (Dirichlet) boundary conditions
* `dbcinds`: indices of the boundary nodes
* `dbcvals`: values of the boundary nodes
* `invinds`: indices of the inner nodes
* `fv`: right hand side of the momentum equation
* `fp`: right hand side of the continuity equation
* `charlen`: characteristic length of the setup
* `nu`: the kinematic viscosity
* `Re`: the Reynolds number
* `odcoo`: dictionary with the coordinates of the domain of \
observation
* `cdcoo`: dictionary with the coordinates of the domain of \
* `ppin` : {int, None}
which dof of `p` is used to pin the pressure, typically \
`-1` for internal flows, and `None` for flows with outflow
control
stokesmatsc : dict
a dictionary of the condensed matrices:
* `M`: the mass matrix of the velocity space,
* `MP`: the mass matrix of the pressure space,
* `A`: the stiffness matrix,
* `JT`: the gradient matrix, and
* `J`: the divergence matrix
* `Jfull`: the uncondensed divergence matrix
and, if `bccontrol=True`, the boundary control matrices that weakly \
impose `Arob*v = Brob*u`, where
* `Arob`: contribution to `A`
* `Brob`: input operator
`if mergerhs`
rhsd : dict
`rhsd_vfrc` and `rhsd_stbc` merged
`else`
rhsd_vfrc : dict
of the dirichlet and pressure fix reduced right hand sides
rhsd_stbc : dict
of the contributions of the boundary data to the rhs:
* `fv`: contribution to momentum equation,
* `fp`: contribution to continuity equation
Examples
--------
femp, stokesmatsc, rhsd_vfrc, rhsd_stbc \
= get_sysmats(problem='drivencavity', N=10, nu=1e-2)
"""
problemdict = dict(drivencavity=drivcav_fems,
cylinderwake=cyl_fems,
# cylinderwake3D=cyl3D_fems,
cylinderwake3D=gen_bccont_fems_3D,
gen_bccont=gen_bccont_fems)
if problem == 'cylinderwake' or problem == 'gen_bccont':
meshparams.update(dict(inflowvel=charvel))
if problem == 'drivencavity':
meshparams = dict(N=meshparams['N'])
if problem == 'cylinder_rot':
problemfem = gen_bccont_fems
meshparams.update(dict(movingwallcntrl=True))
meshparams.update(dict(inflowvel=charvel))
else:
problemfem = problemdict[problem]
femp = problemfem(scheme=scheme, bccontrol=bccontrol, **meshparams)
if onlymesh:
return femp
# setting some parameters
if Re is not None:
nu = charvel*femp['charlen']/Re
else:
Re = charvel*femp['charlen']/nu
if bccontrol:
try:
cbshapefuns = femp['contrbcsshapefuns']
cbclist = femp['contrbcssubdomains']
cbds = None
except KeyError:
cbshapefuns = femp['contrbcsshapefuns']
cbclist = None
cbds = femp['cntrbcsds']
else:
cbclist, cbshapefuns, cbds = None, None, None
try:
outflowds = femp['outflowds']
except KeyError:
outflowds = None
stokesmats = dts.get_stokessysmats(femp['V'], femp['Q'], nu,
cbclist=cbclist, cbds=cbds,
gradvsymmtrc=gradvsymmtrc,
outflowds=outflowds,
cbshapefuns=cbshapefuns,
bccontrol=bccontrol)
rhsd_vf = dts.setget_rhs(femp['V'], femp['Q'],
femp['fv'], femp['fp'], t=0)
# remove the freedom in the pressure if required
if problem == 'cylinderwake':
logging.debug('cylinderwake: pressure need not be pinned')
if ppin is not None:
raise UserWarning('pinning the p will give wrong results')
elif ppin is None:
logging.debug('pressure is not pinned - ' +
'`J` may be singular for internal flow')
elif ppin == -1:
stokesmats['J'] = stokesmats['J'][:-1, :][:, :]
stokesmats['JT'] = stokesmats['JT'][:, :-1][:, :]
rhsd_vf['fp'] = rhsd_vf['fp'][:-1, :]
logging.info('pressure pinned at last dof `-1`')
else:
raise NotImplementedError('Cannot pin `p` other than at `-1`')
# reduce the matrices by resolving the BCs
try:
(stokesmatsc, rhsd_stbc, invinds, _,
_) = dts.condense_sysmatsbybcs(stokesmats, femp['diribcs'])
except KeyError: # can also just provide indices and values
(stokesmatsc, rhsd_stbc, invinds, _, _) = \
dts.condense_sysmatsbybcs(stokesmats, dbcinds=femp['dbcinds'],
dbcvals=femp['dbcvals'])
stokesmatsc.update({'Jfull': stokesmats['J']})
# pressure freedom and dirichlet reduced rhs
rhsd_vfrc = dict(fp=rhsd_vf['fp'], fv=rhsd_vf['fv'][invinds, ])
if bccontrol:
Arob, fvrob = dts.condense_velmatsbybcs(stokesmats['amatrob'],
dbcinds=femp['dbcinds'],
dbcvals=femp['dbcvals'])
if np.linalg.norm(fvrob) > 1e-15:
raise UserWarning('diri and control bc must not intersect')
Brob = stokesmats['bmatrob'][invinds, :]
stokesmatsc.update({'Brob': Brob, 'Arob': Arob})
# add the info on boundary and inner nodes
bcdata = {'invinds': invinds,
'ppin': ppin}
femp.update(bcdata)
femp.update({'nu': nu})
femp.update({'Re': Re})
if mergerhs:
rhsd = dict(fv=rhsd_vfrc['fv']+rhsd_stbc['fv'],
fp=rhsd_vfrc['fp']+rhsd_stbc['fp'])
return femp, stokesmatsc, rhsd
else:
return femp, stokesmatsc, rhsd_vfrc, rhsd_stbc
[docs]def drivcav_fems(N=10, vdgree=2, pdgree=1, scheme=None, bccontrol=None):
"""dictionary for the fem items of the (unit) driven cavity
Parameters
----------
N : int
mesh parameter for the unitsquare (N gives 2*N*N triangles)
vdgree : int, optional
polynomial degree of the velocity basis functions, defaults to 2
pdgree : int, optional
polynomial degree of the pressure basis functions, defaults to 1
scheme : {None, 'CR', 'TH'}
the finite element scheme to be applied, 'CR' for Crouzieux-Raviart,\
'TH' for Taylor-Hood, overrides `pdgree`, `vdgree`, defaults to `None`
bccontrol : boolean, optional
whether to consider boundary control via penalized Robin \
defaults to false. TODO: not implemented yet but we need it here \
for consistency
Returns
-------
femp : a dict
of problem FEM description with the keys:
* `V`: FEM space of the velocity
* `Q`: FEM space of the pressure
* `diribcs`: list of the (Dirichlet) boundary conditions
* `fv`: right hand side of the momentum equation
* `fp`: right hand side of the continuity equation
* `charlen`: characteristic length of the setup
* `odcoo`: dictionary with the coordinates of the domain of \
observation
* `cdcoo`: dictionary with the coordinates of the domain of \
control
"""
mesh = dolfin.UnitSquareMesh(N, N)
if scheme == 'CR':
# print 'we use Crouzieux-Raviart elements !'
V = dolfin.VectorFunctionSpace(mesh, "CR", 1)
Q = dolfin.FunctionSpace(mesh, "DG", 0)
if scheme == 'TH':
# print 'we use Taylor-Hood elements !'
V = dolfin.VectorFunctionSpace(mesh, "CG", 2)
Q = dolfin.FunctionSpace(mesh, "CG", 1)
else:
V = dolfin.VectorFunctionSpace(mesh, "CG", vdgree)
Q = dolfin.FunctionSpace(mesh, "CG", pdgree)
if bccontrol:
raise NotImplementedError()
else:
pass
# Boundaries
def top(x, on_boundary):
return x[1] > 1.0 - dolfin.DOLFIN_EPS
def leftbotright(x, on_boundary):
return (x[0] > 1.0 - dolfin.DOLFIN_EPS
or x[1] < dolfin.DOLFIN_EPS
or x[0] < dolfin.DOLFIN_EPS)
# No-slip boundary condition for velocity
noslip = dolfin.Constant((0.0, 0.0))
bc0 = dolfin.DirichletBC(V, noslip, leftbotright)
# Boundary condition for velocity at the lid
lid = dolfin.Constant(("1", "0.0"))
bc1 = dolfin.DirichletBC(V, lid, top)
# Collect boundary conditions
diribcs = [bc0, bc1]
# rhs of momentum eqn
fv = dolfin.Constant((0.0, 0.0))
# rhs of the continuity eqn
fp = dolfin.Constant(0.0)
dfems = dict(V=V,
Q=Q,
diribcs=diribcs,
fv=fv,
fp=fp,
uspacedep=0,
charlen=1.0)
# domains of observation and control
odcoo = dict(xmin=0.45,
xmax=0.55,
ymin=0.5,
ymax=0.7)
cdcoo = dict(xmin=0.4,
xmax=0.6,
ymin=0.2,
ymax=0.3)
dfems.update(dict(cdcoo=cdcoo, odcoo=odcoo))
return dfems
[docs]def cyl_fems(refinement_level=2, vdgree=2, pdgree=1, scheme=None,
inflowvel=1., bccontrol=False, verbose=False):
"""
dictionary for the fem items for the cylinder wake
Parameters
----------
vdgree : polynomial degree of the velocity basis functions,
defaults to 2
pdgree : polynomial degree of the pressure basis functions,
defaults to 1
scheme : {None, 'CR', 'TH'}
the finite element scheme to be applied, 'CR' for Crouzieux-Raviart,\
'TH' for Taylor-Hood, overrides `pdgree`, `vdgree`, defaults to `None`
bccontrol : boolean, optional
whether to consider boundary control via penalized Robin \
defaults to `False`
Returns
-------
femp : a dictionary with the keys:
* `V`: FEM space of the velocity
* `Q`: FEM space of the pressure
* `diribcs`: list of the (Dirichlet) boundary conditions
* `dirip`: list of the (Dirichlet) boundary conditions \
for the pressure
* `fv`: right hand side of the momentum equation
* `fp`: right hand side of the continuity equation
* `charlen`: characteristic length of the setup
* `odcoo`: dictionary with the coordinates of the \
domain of observation
* `cdcoo`: dictionary with the coordinates of the domain of control
* `uspacedep`: int that specifies in what spatial direction Bu \
changes. The remaining is constant
* `bcsubdoms`: list of subdomains that define the segments where \
the boundary control is applied
Notes
-----
TODO: `inflowvel` as input is there for consistency but not processed
parts of the code were taken from the NSbench collection
https://launchpad.net/nsbench
| __author__ = "Kristian Valen-Sendstad <kvs@simula.no>"
| __date__ = "2009-10-01"
| __copyright__ = "Copyright (C) 2009-2010 " + __author__
| __license__ = "GNU GPL version 3 or any later version"
"""
# Constants related to the geometry
bmarg = 1.e-3 + dolfin.DOLFIN_EPS
xmin = 0.0
xmax = 2.2
ymin = 0.0
ymax = 0.41
xcenter = 0.2
ycenter = 0.2
radius = 0.05
# boundary control at the cylinder
# we define two symmetric wrt x-axis little outlets
# via the radian of the center of the outlets and the extension
centerrad = np.pi/3 # measured from the most downstream point (pi/2 = top)
extensrad = np.pi/6 # radian measure of the extension of the outlets
# bounding boxes for outlet domains
if centerrad + extensrad/2 > np.pi/2 or centerrad - extensrad/2 < 0:
raise NotImplementedError('1st outlet must lie in the 1st quadrant')
b1xmin = xcenter + radius*np.cos(centerrad + extensrad/2)
b1ymax = ycenter + radius*np.sin(centerrad + extensrad/2)
b1xmax = xcenter + radius*np.cos(centerrad - extensrad/2)
b1ymin = ycenter + radius*np.sin(centerrad - extensrad/2)
# symmetry wrt x-axis
b2xmin, b2xmax = b1xmin, b1xmax
b2ymin = ycenter - radius*np.sin(centerrad + extensrad/2)
b2ymax = ycenter - radius*np.sin(centerrad - extensrad/2)
# vectors from the center to the control domain corners
# we need them to define/parametrize the control shape functions
b1base = np.array([[b1xmax - xcenter], [b1ymin - ycenter]])
b2base = np.array([[b2xmin - xcenter], [b2ymin - ycenter]])
# normal vectors of the control domain (considered as a straight line)
centvec = np.array([[xcenter], [ycenter]])
b1tang = np.array([[b1xmax - b1xmin], [b1ymin - b1ymax]])
b2tang = np.array([[b2xmin - b2xmax], [b2ymin - b2ymax]])
rotby90 = np.array([[0, -1.], [1., 0]])
b1normal = rotby90.dot(b1tang) / np.linalg.norm(b1tang)
b2normal = rotby90.dot(b2tang) / np.linalg.norm(b2tang)
if verbose:
print('centvec :', centvec.flatten(), ' b1base', b1base.flatten())
print(b1xmin, b1xmax, b1ymin, b1ymax)
print(b2xmin, b2xmax, b2ymin, b2ymax)
print(b1base, np.linalg.norm(b1base))
print(b1tang)
print(b1normal)
print(b2base, np.linalg.norm(b2base))
print(b2tang)
print(b2normal)
print('diameter of the outlet', radius*2*np.sin(extensrad/2))
print('midpoint of the outlet 1 secant: [{0}, {1}]'.
format(centvec[0]+radius*np.cos(centerrad),
centvec[1]+radius*np.sin(centerrad)))
print('midpoint of the outlet 2 secant: [{0}, {1}]'.
format(centvec[0]+radius*np.cos(centerrad),
centvec[1]-radius*np.sin(centerrad)))
print('angle of midpoint vec 1 and x-axis', np.rad2deg(centerrad))
def insidebbox(x, whichbox=None):
inbbone = (x[0] > b1xmin and x[0] < b1xmax
and x[1] > b1ymin and x[1] < b1ymax)
inbbtwo = (x[0] > b2xmin and x[0] < b2xmax
and x[1] > b2ymin and x[1] < b2ymax)
if whichbox is None:
return inbbone or inbbtwo
if whichbox == 1:
return inbbone
if whichbox == 2:
return inbbtwo
# Inflow boundary
class InflowBoundary(dolfin.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] < xmin + bmarg
# No-slip boundary
class NoslipChannelWalls(dolfin.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (x[1] < ymin + bmarg or x[1] > ymax - bmarg)
class NoslipCylinderSurface(dolfin.SubDomain):
def inside(self, x, on_boundary):
dx = x[0] - xcenter
dy = x[1] - ycenter
r = dolfin.sqrt(dx*dx + dy*dy)
if bccontrol:
notinbbx = not insidebbox(x)
return on_boundary and r < radius + bmarg and notinbbx
else:
return on_boundary and r < radius + bmarg
# Outflow boundary
class OutflowBoundary(dolfin.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] > xmax - bmarg
# # meshes are available from https://launchpad.net/nsbench
# Load mesh
if refinement_level > 9:
raise RuntimeError("No mesh available for refinement level {0}".
format(refinement_level))
try:
mesh = dolfin.Mesh("mesh/cylinder_%d.xml.gz" % refinement_level)
except RuntimeError:
mesh = dolfin.Mesh("mesh/cylinder_%d.xml" % refinement_level)
# scheme = 'CR'
if scheme == 'CR':
# print 'we use Crouzieux-Raviart elements !'
V = dolfin.VectorFunctionSpace(mesh, "CR", 1)
Q = dolfin.FunctionSpace(mesh, "DG", 0)
else:
V = dolfin.VectorFunctionSpace(mesh, "CG", vdgree)
Q = dolfin.FunctionSpace(mesh, "CG", pdgree)
if bccontrol:
def _csf(s, nvec):
return ((1. - 0.5*(1 + np.sin(s*2*np.pi + 0.5*np.pi)))*nvec[0],
(1. - 0.5*(1 + np.sin(s*2*np.pi + 0.5*np.pi)))*nvec[1])
class ContBoundaryOne(dolfin.SubDomain):
def inside(self, x, on_boundary):
dx = x[0] - xcenter
dy = x[1] - ycenter
r = dolfin.sqrt(dx*dx + dy*dy)
inbbx = insidebbox(x, whichbox=1)
return on_boundary and r < radius + bmarg and inbbx
def cont_shape_one(element=None):
class ContShapeOne(dolfin.UserExpression):
def eval(self, value, x):
xvec = x - centvec.flatten()
aang = np.arccos(np.dot(xvec, b1base)
/ (np.linalg.norm(xvec)
* np.linalg.norm(b1base)))
s = aang/extensrad
vls = _csf(s, b1normal)
value[0], value[1] = vls[0], vls[1]
if verbose:
dx = x[0] - xcenter
dy = x[1] - ycenter
r = dolfin.sqrt(dx*dx + dy*dy)
print(x - centvec.flatten(), ': s=', s, ': r=', r,
':', np.linalg.norm(np.array(vls)))
def value_shape(self):
return (2,)
return ContShapeOne(element=element)
class ContBoundaryTwo(dolfin.SubDomain):
def inside(self, x, on_boundary):
dx = x[0] - xcenter
dy = x[1] - ycenter
r = dolfin.sqrt(dx*dx + dy*dy)
inbbx = insidebbox(x, whichbox=2)
return on_boundary and r < radius + bmarg and inbbx
def cont_shape_two(element=None):
class ContShapeTwo(dolfin.UserExpression):
def eval(self, value, x):
xvec = x - centvec.flatten()
aang = np.arccos(np.dot(xvec, b2base)
/ (np.linalg.norm(xvec)
* np.linalg.norm(b2base)))
s = aang/extensrad
vls = _csf(s, b2normal)
value[0], value[1] = vls[0], vls[1]
if verbose:
dx = x[0] - xcenter
dy = x[1] - ycenter
r = dolfin.sqrt(dx*dx + dy*dy)
print(x - centvec.flatten(), ': s=', s, ': r=', r,
':', np.linalg.norm(np.array(vls)))
def value_shape(self):
return (2,)
return ContShapeTwo(element=element)
bcsubdoms = [ContBoundaryOne, ContBoundaryTwo]
bcshapefuns = [cont_shape_one(element=V.ufl_element()),
cont_shape_two(element=V.ufl_element())]
else:
bcsubdoms = [None, None]
bcshapefuns = [None, None]
# dolfin.plot(mesh)
# dolfin.interactive(True)
# Create right-hand side function
fv = dolfin.Constant((0, 0))
fp = dolfin.Constant(0)
def initial_conditions(self, V, Q):
u0 = dolfin.Constant((0, 0))
p0 = dolfin.Constant(0)
return u0, p0
# Create inflow boundary condition
g0 = dolfin.Expression(('4*(x[1]*(ymax-x[1]))/(ymax*ymax)', '0.0'),
ymax=ymax, element=V.ufl_element())
bc0 = dolfin.DirichletBC(V, g0, InflowBoundary())
# Create no-slip boundary condition
g1 = dolfin.Constant((0, 0))
bc1 = dolfin.DirichletBC(V, g1, NoslipChannelWalls())
# Create no-slip at cylinder surface
bc1cyl = dolfin.DirichletBC(V, g1, NoslipCylinderSurface())
# Create outflow boundary condition for pressure
g2 = dolfin.Constant(0)
bc2 = dolfin.DirichletBC(Q, g2, OutflowBoundary())
# Collect boundary conditions
bcu = [bc0, bc1, bc1cyl]
bcp = [bc2]
dbcinds, dbcvals = [], []
for bc in bcu:
bcdict = bc.get_boundary_values()
dbcvals.extend(list(bcdict.values()))
dbcinds.extend(list(bcdict.keys()))
cylfems = dict(V=V,
Q=Q,
diribcs=bcu,
dbcinds=dbcinds,
dbcvals=dbcvals,
dirip=bcp,
contrbcssubdomains=bcsubdoms,
contrbcsshapefuns=bcshapefuns,
fv=fv,
fp=fp,
uspacedep=0,
charlen=0.1,
mesh=mesh)
# domains of observation and control
odcoo = dict(xmin=0.6,
xmax=0.7,
ymin=0.15,
ymax=0.25)
cdcoo = dict(xmin=0.27,
xmax=0.32,
ymin=0.15,
ymax=0.25)
cylfems.update(dict(cdcoo=cdcoo, odcoo=odcoo))
return cylfems
[docs]def cyl3D_fems(refinement_level=2, scheme='TH',
strtobcsobs='',
strtomeshfile='', strtophysicalregions='',
bccontrol=False, verbose=False):
"""
dictionary for the fem items for the 3D cylinder wake
which is
* the 2D setup extruded in z-direction
* with symmetry BCs at the z-walls
Parameters
----------
scheme : {None, 'CR', 'TH'}
the finite element scheme to be applied, 'CR' for Crouzieux-Raviart,\
'TH' for Taylor-Hood, overrides `pdgree`, `vdgree`, defaults to `None`
bccontrol : boolean, optional
whether to consider boundary control via penalized Robin \
defaults to `False`
Returns
-------
femp : a dictionary with the keys:
* `V`: FEM space of the velocity
* `Q`: FEM space of the pressure
* `diribcs`: list of the (Dirichlet) boundary conditions
* `dirip`: list of the (Dirichlet) boundary conditions \
for the pressure
* `fv`: right hand side of the momentum equation
* `fp`: right hand side of the continuity equation
* `charlen`: characteristic length of the setup
* `odcoo`: dictionary with the coordinates of the \
domain of observation
* `cdcoo`: dictionary with the coordinates of the domain of control
* `uspacedep`: int that specifies in what spatial direction Bu \
changes. The remaining is constant
* `bcsubdoms`: list of subdomains that define the segments where \
the boundary control is applied
Notes
-----
parts of the code were taken from the NSbench collection
https://launchpad.net/nsbench
| __author__ = "Kristian Valen-Sendstad <kvs@simula.no>"
| __date__ = "2009-10-01"
| __copyright__ = "Copyright (C) 2009-2010 " + __author__
| __license__ = "GNU GPL version 3 or any later version"
"""
# Constants related to the geometry
# xmin = 0.0
# xmax = 8.0
# ymin = 0.0
ymax = 1.5
# zmin = 0.0
# zmax = 1.0
# xcenter = 2.0
# ycenter = ymax/2
# radius = 0.15
# Load mesh
if not strtomeshfile == '':
mesh = dolfin.Mesh(strtomeshfile)
meshfile = strtophysicalregions
else:
mesh = dolfin.\
Mesh("mesh/3d-cyl/karman3D_lvl{0}.xml.gz".format(refinement_level))
meshfile = 'mesh/3d-cyl/karman3D_lvl{0}_facet_region.xml.gz'.\
format(refinement_level)
# scheme = 'CR'
if scheme == 'CR':
# print 'we use Crouzieux-Raviart elements !'
V = dolfin.VectorFunctionSpace(mesh, "CR", 1)
Q = dolfin.FunctionSpace(mesh, "DG", 0)
elif scheme == 'TH':
V = dolfin.VectorFunctionSpace(mesh, "CG", 2)
Q = dolfin.FunctionSpace(mesh, "CG", 1)
# get the boundaries from the gmesh file
boundaries = dolfin.\
MeshFunction('size_t', mesh, meshfile)
# # Create inflow boundary condition
# if no-slip at front and back
# gin = dolfin.Expression(('6*(x[1]*(ymax-x[1]))/(ymax*ymax) * ' +
# '6*(x[2]*(zmax-x[2]))/(zmax*zmax)',
# '0.0', '0.0'),
# ymax=ymax, zmax=zmax, element=V.ufl_element())
gin = dolfin.Expression(('6*(x[1]*(ymax-x[1]))/(ymax*ymax)', '0.0', '0.0'),
ymax=ymax, element=V.ufl_element())
bcin = dolfin.DirichletBC(V, gin, boundaries, 1)
# ## Create no-slip boundary condition
gzero = dolfin.Constant((0, 0, 0))
# channel walls
bcbt = dolfin.DirichletBC(V, gzero, boundaries, 2) # bottom
bctp = dolfin.DirichletBC(V, gzero, boundaries, 6) # top
# yes-slip at front and back -- only z-component set to zero
gscalzero = dolfin.Constant(0)
bcbc = dolfin.DirichletBC(V.sub(2), gscalzero, boundaries, 5) # back
bcfr = dolfin.DirichletBC(V.sub(2), gscalzero, boundaries, 4) # front
# Create no-slip at cylinder surface
bccuc = dolfin.DirichletBC(V, gzero, boundaries, 9) # uncontrolled
bccco = dolfin.DirichletBC(V, gzero, boundaries, 7) # ctrl upper
bccct = dolfin.DirichletBC(V, gzero, boundaries, 8) # ctrl lower
# Create outflow boundary condition for pressure
g2 = dolfin.Constant(0)
bc2 = dolfin.DirichletBC(Q, g2, boundaries, 3)
# Collect boundary conditions
bcu = [bcin, bcbt, bcfr, bcbc, bctp, bccuc, bccco, bccct]
bcp = [bc2]
# Create right-hand side function
fv = dolfin.Constant((0, 0, 0))
fp = dolfin.Constant(0)
def initial_conditions(self, V, Q):
u0 = dolfin.Constant((0, 0, 0))
p0 = dolfin.Constant(0)
return u0, p0
cylfems = dict(V=V,
Q=Q,
diribcs=bcu,
dirip=bcp,
# contrbcssubdomains=bcsubdoms,
# contrbcsshapefuns=bcshapefuns,
fv=fv,
fp=fp,
uspacedep=0,
charlen=0.3,
mesh=mesh)
return cylfems
[docs]def gen_bccont_fems(scheme='TH', bccontrol=True, verbose=False,
strtomeshfile='', strtophysicalregions='',
inflowvel=1., inflowprofile='parabola',
movingwallcntrl=False,
strtobcsobs=''):
"""
dictionary for the fem items for a general 2D flow setup
with
* inflow/outflow
* boundary control
Parameters
----------
scheme : {None, 'CR', 'TH'}
the finite element scheme to be applied, 'CR' for Crouzieux-Raviart,\
'TH' for Taylor-Hood, overrides `pdgree`, `vdgree`, defaults to `None`
bccontrol : boolean, optional
whether to consider boundary control via penalized Robin \
defaults to `True`
movingwallcntrl : boolean, optional
whether control is via moving boundaries
Returns
-------
femp : a dictionary with the keys:
* `V`: FEM space of the velocity
* `Q`: FEM space of the pressure
* `diribcs`: list of the (Dirichlet) boundary conditions
* `dbcsinds`: list vortex indices with (Dirichlet) boundary conditions
* `dbcsvals`: list of values of the (Dirichlet) boundary conditions
* `dirip`: list of the (Dirichlet) boundary conditions \
for the pressure
* `fv`: right hand side of the momentum equation
* `fp`: right hand side of the continuity equation
* `charlen`: characteristic length of the setup
* `odcoo`: dictionary with the coordinates of the \
domain of observation
"""
# Load mesh
logging.info('mesh: ' + strtomeshfile)
mesh = dolfin.Mesh(strtomeshfile)
if scheme == 'CR':
V = dolfin.VectorFunctionSpace(mesh, "CR", 1)
Q = dolfin.FunctionSpace(mesh, "DG", 0)
elif scheme == 'TH':
V = dolfin.VectorFunctionSpace(mesh, "CG", 2)
Q = dolfin.FunctionSpace(mesh, "CG", 1)
else:
raise UserWarning('Specify the FEM scheme `{TH, CR}` explicitly SVP')
boundaries = dolfin.MeshFunction('size_t', mesh, strtophysicalregions)
with open(strtobcsobs) as f:
cntbcsdata = json.load(f)
inflowgeodata = cntbcsdata['inflow']
inflwpe = inflowgeodata['physical entity']
logging.info('mesh: physical entity {0} -- inflow'.format(inflwpe))
inflwin = np.array(inflowgeodata['inward normal'])
inflwxi = np.array(inflowgeodata['xone'])
inflwxii = np.array(inflowgeodata['xtwo'])
leninflwb = np.linalg.norm(inflwxi-inflwxii)
if inflowprofile == 'block':
inflwprfl = dolfin.\
Expression(('cv*no', 'cv*nt'), cv=inflowvel,
no=inflwin[0], nt=inflwin[1],
element=V.ufl_element())
elif inflowprofile == 'parabola':
inflwprfl = InflowParabola(degree=2, lenb=leninflwb, xone=inflwxi,
normalvec=inflwin, inflowvel=inflowvel)
bcin = dolfin.DirichletBC(V, inflwprfl, boundaries, inflwpe)
diribcu = [bcin]
# ## THE WALLS
wallspel = cntbcsdata['walls']['physical entity']
gzero = dolfin.Constant((0, 0))
for wpe in wallspel:
diribcu.append(dolfin.DirichletBC(V, gzero, boundaries, wpe))
bcdict = diribcu[-1].get_boundary_values()
logging.info('mesh: physical entity {0} -- wall'.format(wpe))
if not bccontrol: # treat the control boundaries as walls
try:
for cntbc in cntbcsdata['controlbcs']:
diribcu.append(dolfin.DirichletBC(V, gzero, boundaries,
cntbc['physical entity']))
logging.info('mesh: physical entity {0} -- wall'.
format(cntbc['physical entity']))
except KeyError:
pass # no control boundaries
mvwdbcs = []
mvwtvs = []
try:
for cntbc in cntbcsdata['moving walls']:
if cntbc['type'] == 'circle':
center = np.array(cntbc['geometry']['center'])
radius = cntbc['geometry']['radius']
omega = 1. if movingwallcntrl else 0.
rotcyl = RotatingCircle(degree=2, radius=radius,
center=center, omega=omega)
else:
raise NotImplementedError()
mvwdbcs.append(dolfin.DirichletBC(V, rotcyl, boundaries,
cntbc['physical entity']))
logging.info('mesh: physical entity {0} -- moving wall'.
format(cntbc['physical entity']))
except KeyError:
pass # no moving walls defined
if not movingwallcntrl and len(mvwdbcs) > 0:
diribcu.extend(mvwdbcs) # add the moving walls to the diri bcs
mvwdbcs = []
logging.info('mesh: physical entities -- moving walls --> walls')
# Create outflow boundary condition for pressure
# TODO XXX why zero pressure?? is this do-nothing???
outflwpe = cntbcsdata['outflow']['physical entity']
logging.info('mesh: physical entity {0} -- outflow'.format(outflwpe))
g2 = dolfin.Constant(0)
bc2 = dolfin.DirichletBC(Q, g2, boundaries, outflwpe)
# Collect boundary conditions
bcp = [bc2]
# Create right-hand side function
fv = dolfin.Constant((0, 0))
fp = dolfin.Constant(0)
def initial_conditions(self, V, Q):
u0 = dolfin.Constant((0, 0))
p0 = dolfin.Constant(0)
return u0, p0
dbcinds, dbcvals = [], []
for bc in diribcu:
bcdict = bc.get_boundary_values()
dbcvals.extend(list(bcdict.values()))
dbcinds.extend(list(bcdict.keys()))
mvwbcinds, mvwbcvals = [], []
for bc in mvwdbcs:
bcdict = bc.get_boundary_values()
mvwbcvals.extend(list(bcdict.values()))
mvwbcinds.extend(list(bcdict.keys()))
# ## Control boundaries
bcpes, bcshapefuns, bcds = [], [], []
if bccontrol:
for cbc in cntbcsdata['controlbcs']:
if cbc['type'] == 'inlet':
cxi, cxii = np.array(cbc['xone']), np.array(cbc['xtwo'])
csf = _get_cont_shape_fun2D(xi=cxi, xii=cxii,
element=V.ufl_element())
elif cbc['type'] == 'rotating circle':
csf = RotatingCircle(center=cbc['center'],
radius=cbc['radius'])
cpe = cbc['physical entity']
logging.info('mesh: physical entity {0} -- boundary control ({1})'.
format(cpe, cbc['type']))
bcshapefuns.append(csf)
bcpes.append(cpe)
bcds.append(dolfin.Measure("ds", subdomain_data=boundaries)(cpe))
# ## Lift Drag Computation
try:
ldsurfpe = cntbcsdata['lift drag surface']['physical entity']
liftdragds = dolfin.Measure("ds", subdomain_data=boundaries)(ldsurfpe)
bclds = dolfin.DirichletBC(V, gzero, boundaries, ldsurfpe)
logging.info(f'mesh: physical entity {ldsurfpe} -- lift/drag surf')
bcldsdict = bclds.get_boundary_values()
ldsbcinds = list(bcldsdict.keys())
except KeyError:
liftdragds = None # no domain specified for lift/drag
ldsbcinds = None
try:
outflwpe = cntbcsdata['outflow']['physical entity']
outflowds = dolfin.Measure("ds", subdomain_data=boundaries)(outflwpe)
except KeyError:
outflowds = None # no domain specified for outflow
try:
odcoo = cntbcsdata['observation-domain-coordinates']
except KeyError:
odcoo = None
gbcfems = dict(V=V,
Q=Q,
dbcinds=dbcinds,
dbcvals=dbcvals,
mvwbcinds=mvwbcinds,
mvwbcvals=mvwbcvals,
mvwtvs=mvwtvs,
dirip=bcp,
outflowds=outflowds,
# contrbcssubdomains=bcsubdoms,
liftdragds=liftdragds,
ldsbcinds=ldsbcinds,
contrbcmeshfunc=boundaries,
contrbcspes=bcpes,
contrbcsshapefuns=bcshapefuns,
cntrbcsds=bcds,
odcoo=odcoo,
fv=fv,
fp=fp,
charlen=cntbcsdata['characteristic length'],
mesh=mesh)
return gbcfems
def _get_cont_shape_fun2D(xi=None, xii=None, element=None, shape='parabola'):
lencb = np.linalg.norm(xi-xii)
cbt = 1./lencb*(xii-xi) # the normalized vector pointing x1 -> x2
cbn = np.array([cbt[1], -cbt[0]]).reshape((2, 1)) # rotate by pi/2
class GenContShape(dolfin.UserExpression):
def __init__(self, degree=2, element=None):
self.degree = degree
self.element = element
super().__init__()
def eval(self, value, x):
curs = np.linalg.norm(x - xi)/lencb
logging.debug(f'x={x}: s={curs}')
curvel = 6*curs*(1-curs)*cbn
value[0], value[1] = curvel[0], curvel[1]
def value_shape(self):
return (2,)
return GenContShape(element=element)
class InflowParabola(dolfin.UserExpression):
'''Create inflow boundary condition
a parabola g with `int g(s)ds = s1-s0 == int 1 ds`
if on [0,1]: `g(s) = s*(1-s)*4*3/2`
if on [s0,s1]: `g(s) = ((s-s0)/(s1-s0))*(1-(s-s0)/(s1-s0))*6`
since then `g((s0+s1)/2)=3/2` and `g(s0)=0=g(s1)`'''
def __init__(self, degree=2, lenb=None, xone=None,
inflowvel=1., normalvec=None):
self.degree = degree
self.lenb = lenb
self.xone = xone
self.normalvec = normalvec
self.inflowvel = inflowvel
try:
super().__init__()
except RuntimeError():
pass # had trouble with this call to __init__ in 'dolfin:2017.2.0'
def eval(self, value, x):
curs = np.linalg.norm(x - self.xone)/self.lenb
# print(x, curs)
curvel = self.inflowvel*6*curs*(1-curs)*self.normalvec
value[0], value[1] = curvel[0], curvel[1]
def value_shape(self):
return (2,)
class InflowParabola3D(dolfin.UserExpression):
'''Create inflow boundary condition
a parabola g with `int g(s)ds = s1-s0 == int 1 ds`
if on [0,1]: `g(s) = s*(1-s)*4*3/2`
if on [s0,s1]: `g(s) = ((s-s0)/(s1-s0))*(1-(s-s0)/(s1-s0))*6`
since then `g((s0+s1)/2)=3/2` and `g(s0)=0=g(s1)`'''
def __init__(self, degree=2, xone=None, xtwo=None, xfour=None,
inflowvel=1., normalvec=None):
self.degree = degree
self.xone = xone
self.normalvec = normalvec
self.inflowvel = inflowvel
self.xvec = xtwo-xone
self.yvec = xfour-xone
self.lenxsqrd = np.inner(self.xvec, self.xvec)
self.lenysqrd = np.inner(self.yvec, self.yvec)
try:
super().__init__()
except RuntimeError():
pass # had trouble with this call to __init__ in 'dolfin:2017.2.0'
def eval(self, value, x):
xclean = x - self.xone
cursx = np.inner(xclean, self.xvec)/self.lenxsqrd
cursy = np.inner(xclean, self.yvec)/self.lenysqrd
# print(x, cursx, cursy)
# gin = dolfin.Expression(('6*(x[1]*(ymax-x[1]))/(ymax*ymax)',
# '0.0', '0.0'),
# ymax=ymax, element=V.ufl_element())
curvel = self.inflowvel*36*cursx*(1-cursx) * \
cursy*(1-cursy)*self.normalvec
value[0], value[1], value[2] = curvel[0], curvel[1], curvel[2]
def value_shape(self):
return (3,)
class RotatingCircle(dolfin.UserExpression):
'''Create the boundary condition of a rotating circle
returns the angular velocity at the circle boundary
'''
def __init__(self, degree=2, radius=None, center=None,
omega=1.):
self.degree = degree
self.radius = radius
self.xcenter = center
self.anglevel = radius*omega
print('Rotating cylinder: omega set to {0}'.format(omega))
super().__init__()
def eval(self, value, x):
curn = 1./self.radius*(x - self.xcenter)
# print(np.linalg.norm(curn))
value[0], value[1] = -self.anglevel*curn[1], self.anglevel*curn[0]
def value_shape(self):
return (2,)
class LiftDragSurfForce():
def __init__(self, V=None, nu=None, ldds=None, gradvsymmtrc=True,
outflowds=None, phione=None, phitwo=None):
self.mesh = V.mesh()
self.n = dolfin.FacetNormal(self.mesh)
# self.Id = dolfin.Identity(self.mesh.geometry().dim())
# self.ldds = ldds
self.nu = nu
self.A = dolfin.as_matrix([[0., 1.],
[-1., 0.]])
pickx = dolfin.as_matrix([[1., 0.], [0., 0.]])
picky = dolfin.as_matrix([[0., 0.], [0., 1.]])
self.phionex = pickx*phione
self.phioney = picky*phione
self.phitwo = phitwo
self.gradvsymmtrc = gradvsymmtrc
if gradvsymmtrc:
def epsilon(u):
return 0.5*(dolfin.nabla_grad(u) + dolfin.nabla_grad(u).T)
self.outflowds = outflowds
else:
def epsilon(u):
return dolfin.grad(u)
self.epsilon = epsilon
def evaliftdragforce(self, u=None, p=None):
# ## The direct way
# ux = u.sub(0)
# uy = u.sub(1)
# D = (2*self.nu*inner(self.epsilon(ux), dolfin.grad(self.phitwo))
# # inner(self.nu*grad(ux), grad(self.phione))
# + inner(u, grad(ux))*self.phione
# - p*self.phione.dx(0))*dolfin.dx
# L = (inner(self.nu*grad(uy), grad(self.phione))
# + inner(u, grad(uy))*self.phione
# - p*self.phione.dx(1))*dolfin.dx
# T = -p*self.Id + 2.0*self.nu*self.epsilon(u)
# force = dolfin.dot(T, self.n)
# D = force[0]*self.ldds
# L = force[1]*self.ldds
# drag = dolfin.assemble(L)
# lift = dolfin.assemble(D)
inner, dx = dolfin.inner, dolfin.dx
epsi, grad, div = self.epsilon, dolfin.grad, dolfin.div
# ## the Babuska/Millner trick
pox = self.phionex
poy = self.phioney
# convection
drgo = inner(dolfin.dot(u, dolfin.nabla_grad(u)), pox)*dx
# outflow correction for symmetrized gradient
if self.gradvsymmtrc:
dofcor = (self.nu*inner(grad(u).T*self.n, pox))*self.outflowds
lofcor = (self.nu*inner(grad(u).T*self.n, poy))*self.outflowds
else:
dofcor, lofcor = 0, 0
# diffusion
drgt = (2*self.nu*inner(epsi(u), grad(self.phionex)))*dx
# pressure
drgd = -p*div(self.phionex)*dx
drag = dolfin.assemble(drgt+drgd+drgo-dofcor)
lfto = inner(dolfin.dot(u, dolfin.nabla_grad(u)), poy)*dx
lftt = 2*self.nu*inner(self.epsilon(u), dolfin.grad(poy))*dx
lftd = (-p*dolfin.div(poy))*dx
lift = dolfin.assemble(lfto+lftt+lftd-lofcor)
return lift, drag
def evatorqueSphere2D(self, u=None, p=None):
inner, dx, grad = dolfin.inner, dolfin.dx, dolfin.grad
donto = inner(dolfin.dot(u, dolfin.nabla_grad(u)), self.phitwo)*dx
dontt = 2*self.nu*inner(self.epsilon(u), dolfin.grad(self.phitwo))*dx
dontd = (-p*dolfin.div(self.phitwo))*dx
ofcor = (self.nu*inner(grad(u).T*self.n, self.phitwo))*self.outflowds
tconv = dolfin.assemble(donto)
tdiff = dolfin.assemble(dontt)
tpres = dolfin.assemble(dontd)
tofcr = dolfin.assemble(ofcor)
trqthr = tconv + tdiff + tpres - tofcr
return trqthr
[docs]def gen_bccont_fems_3D(scheme='TH', bccontrol=True, verbose=False,
strtomeshfile='', strtophysicalregions='',
inflowvel=1., inflowprofile='parabola',
movingwallcntrl=False,
strtobcsobs=''):
"""
dictionary for the fem items for a general 3D flow setup
with
* inflow/outflow
* boundary control
Parameters
----------
scheme : {None, 'CR', 'TH'}
the finite element scheme to be applied, 'CR' for Crouzieux-Raviart,\
'TH' for Taylor-Hood, overrides `pdgree`, `vdgree`, defaults to `None`
bccontrol : boolean, optional
whether to consider boundary control via penalized Robin \
defaults to `True`
movingwallcntrl : boolean, optional
whether control is via moving boundaries
Returns
-------
femp : a dictionary with the keys:
* `V`: FEM space of the velocity
* `Q`: FEM space of the pressure
* `diribcs`: list of the (Dirichlet) boundary conditions
* `dbcsinds`: list vortex indices with (Dirichlet) boundary conditions
* `dbcsvals`: list of values of the (Dirichlet) boundary conditions
* `dirip`: list of the (Dirichlet) boundary conditions \
for the pressure
* `fv`: right hand side of the momentum equation
* `fp`: right hand side of the continuity equation
* `charlen`: characteristic length of the setup
* `odcoo`: dictionary with the coordinates of the \
domain of observation
"""
# Load mesh
mesh = dolfin.Mesh(strtomeshfile)
if scheme == 'CR':
V = dolfin.VectorFunctionSpace(mesh, "CR", 1)
Q = dolfin.FunctionSpace(mesh, "DG", 0)
elif scheme == 'TH':
V = dolfin.VectorFunctionSpace(mesh, "CG", 2)
Q = dolfin.FunctionSpace(mesh, "CG", 1)
boundaries = dolfin.MeshFunction('size_t', mesh, strtophysicalregions)
with open(strtobcsobs) as f:
cntbcsdata = json.load(f)
inflowgeodata = cntbcsdata['inflow']
inflwpe = inflowgeodata['physical entity']
inflwin = np.array(inflowgeodata['inward normal'])
inflwxi = np.array(inflowgeodata['xone'])
inflwxii = np.array(inflowgeodata['xtwo'])
# inflwxiii = np.array(inflowgeodata['xthree'])
inflwxiv = np.array(inflowgeodata['xfour'])
if inflowprofile == 'block':
raise NotImplementedError()
inflwprfl = dolfin.\
Expression(('cv*no', 'cv*nt'), cv=inflowvel,
no=inflwin[0], nt=inflwin[1],
element=V.ufl_element())
elif inflowprofile == 'parabola':
inflwprfl = InflowParabola3D(degree=2, xone=inflwxi, xtwo=inflwxii,
xfour=inflwxiv,
normalvec=inflwin, inflowvel=inflowvel)
bcin = dolfin.DirichletBC(V, inflwprfl, boundaries, inflwpe)
diribcu = [bcin]
# ## THE WALLS
wallspel = cntbcsdata['walls']['physical entity']
gzero = dolfin.Constant((0, 0, 0))
for wpe in wallspel:
diribcu.append(dolfin.DirichletBC(V, gzero, boundaries, wpe))
# bcdict = diribcu[-1].get_boundary_values()
if not bccontrol: # treat the control boundaries as walls
try:
for cntbc in cntbcsdata['controlbcs']:
diribcu.append(dolfin.DirichletBC(V, gzero, boundaries,
cntbc['physical entity']))
except KeyError:
pass # no control boundaries
# yes-slip walls
try:
slipwallspel = cntbcsdata['slipwalls']['physical entity']
slipwallsnvs = cntbcsdata['slipwalls']['inward normals']
gscalzero = dolfin.Constant(0)
for kk, swpe in enumerate(slipwallspel):
cinwnrml = np.array(slipwallsnvs[kk])
if np.abs(np.inner(cinwnrml, np.array([0, 0, 1.]))) == 1:
cbcsw = dolfin.DirichletBC(V.sub(2), gscalzero,
boundaries, swpe)
diribcu.append(cbcsw)
else:
raise NotImplementedError()
except KeyError:
pass # no control boundaries
mvwdbcs = []
mvwtvs = []
try:
for cntbc in cntbcsdata['moving walls']:
raise NotImplementedError()
center = np.array(cntbc['geometry']['center'])
radius = cntbc['geometry']['radius']
if cntbc['type'] == 'circle':
omega = 1. if movingwallcntrl else 0.
rotcyl = RotatingCircle(degree=2, radius=radius,
center=center, omega=omega)
else:
raise NotImplementedError()
mvwdbcs.append(dolfin.DirichletBC(V, rotcyl, boundaries,
cntbc['physical entity']))
except KeyError:
pass # no moving walls defined
if not movingwallcntrl:
diribcu.extend(mvwdbcs) # add the moving walls to the diri bcs
mvwdbcs = []
# Create outflow boundary condition for pressure
# TODO XXX why zero pressure?? is this do-nothing???
# outflwpe = cntbcsdata['outflow']['physical entity']
# g2 = dolfin.Constant(0)
# bc2 = dolfin.DirichletBC(Q, g2, boundaries, outflwpe)
# Collect boundary conditions
# bcp = [bc2]
# Create right-hand side function
fv = dolfin.Constant((0, 0, 0))
fp = dolfin.Constant(0)
def initial_conditions(self, V, Q):
u0 = dolfin.Constant((0, 0, 0))
p0 = dolfin.Constant(0)
return u0, p0
dbcinds, dbcvals = [], []
for bc in diribcu:
bcdict = bc.get_boundary_values()
dbcvals.extend(list(bcdict.values()))
dbcinds.extend(list(bcdict.keys()))
mvwbcinds, mvwbcvals = [], []
for bc in mvwdbcs:
bcdict = bc.get_boundary_values()
mvwbcvals.extend(list(bcdict.values()))
mvwbcinds.extend(list(bcdict.keys()))
# ## Control boundaries
bcpes, bcshapefuns, bcds = [], [], []
if bccontrol:
raise NotImplementedError()
for cbc in cntbcsdata['controlbcs']:
cpe = cbc['physical entity']
cxi, cxii = np.array(cbc['xone']), np.array(cbc['xtwo'])
csf = _get_cont_shape_fun2D(xi=cxi, xii=cxii,
element=V.ufl_element())
bcshapefuns.append(csf)
bcpes.append(cpe)
bcds.append(dolfin.Measure("ds", subdomain_data=boundaries)(cpe))
# ## Lift Drag Computation
try:
ldsurfpe = cntbcsdata['lift drag surface']['physical entity']
raise NotImplementedError()
liftdragds = dolfin.Measure("ds", subdomain_data=boundaries)(ldsurfpe)
bclds = dolfin.DirichletBC(V, gzero, boundaries, ldsurfpe)
bcldsdict = bclds.get_boundary_values()
ldsbcinds = list(bcldsdict.keys())
except KeyError:
liftdragds = None # no domain specified for lift/drag
ldsbcinds = None
try:
outflwpe = cntbcsdata['outflow']['physical entity']
outflowds = dolfin.Measure("ds", subdomain_data=boundaries)(outflwpe)
except KeyError:
outflowds = None # no domain specified for outflow
try:
odcoo = cntbcsdata['observation-domain-coordinates']
raise NotImplementedError()
except KeyError:
odcoo = None
gbcfems = dict(V=V,
Q=Q,
dbcinds=dbcinds,
dbcvals=dbcvals,
mvwbcinds=mvwbcinds,
mvwbcvals=mvwbcvals,
mvwtvs=mvwtvs,
# dirip=bcp,
outflowds=outflowds,
# contrbcssubdomains=bcsubdoms,
liftdragds=liftdragds,
ldsbcinds=ldsbcinds,
contrbcmeshfunc=boundaries,
contrbcspes=bcpes,
contrbcsshapefuns=bcshapefuns,
cntrbcsds=bcds,
odcoo=odcoo,
fv=fv,
fp=fp,
charlen=cntbcsdata['characteristic length'],
mesh=mesh)
return gbcfems
def get_bcinds(mesh=None, V=None, prfile=None, pelist=[]):
boundaries = dolfin.MeshFunction('size_t', mesh, prfile)
gzero = dolfin.Constant((0, 0))
diribcinds = []
for wpe in pelist:
cdbc = dolfin.DirichletBC(V, gzero, boundaries, wpe)
bcdict = cdbc.get_boundary_values()
diribcinds.extend(list(bcdict.keys()))
return diribcinds