Skip to content

Ensuring consistency in FEniCS scripts for FSI tutorial cases #142

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Jan 8, 2021
Merged

Ensuring consistency in FEniCS scripts for FSI tutorial cases #142

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8c149d9
Making SU2- and OpenFOAM- case of perp-flap for FEniCS consistent
IshaanDesai Dec 31, 2020
fba42f6
Reverting wrong change in Allclean script
IshaanDesai Dec 31, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 19 additions & 21 deletions FSI/flap_perp_2D/OpenFOAM-FEniCS/Solid/perp-flap.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,32 +1,25 @@
# Import required libs
from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \
TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, SubDomain, \
Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system, MPI, MeshFunction
import dolfin

TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, Identity, inner, dx, ds, \
sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system, MPI, MeshFunction
from ufl import nabla_div
import numpy as np
import matplotlib.pyplot as plt
from fenicsprecice import Adapter
from enum import Enum


class clampedBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
if on_boundary and abs(x[1]) < tol:
return True
else:
return False
# define the two kinds of boundary: clamped and coupling Neumann Boundary
def clamped_boundary(x, on_boundary):
return on_boundary and abs(x[1]) < tol


def neumann_boundary(x, on_boundary):
"""
determines whether a node is on the coupling boundary

class neumannBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
if on_boundary and ((abs(x[1] - 1) < tol) or abs(abs(x[0]) - W / 2) < tol):
return True
else:
return False
"""
return on_boundary and ((abs(x[1] - 1) < tol) or abs(abs(x[0]) - W / 2) < tol)


# Geometry and material properties
Expand All @@ -51,6 +44,9 @@ def inside(self, x, on_boundary):
# create Function Space
V = VectorFunctionSpace(mesh, 'P', 2)

# BCs
tol = 1E-14

# Trial and Test Functions
du = TrialFunction(V)
v = TestFunction(V)
Expand All @@ -66,19 +62,21 @@ def inside(self, x, on_boundary):
f_N_function = interpolate(Expression(("1", "0"), degree=1), V)
u_function = interpolate(Expression(("0", "0"), degree=1), V)

coupling_boundary = AutoSubDomain(neumann_boundary)
fixed_boundary = AutoSubDomain(clamped_boundary)

precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json")

# Initialize the coupling interface
# Function space V is passed twice as both read and write functions are defined using the same space
precice_dt = precice.initialize(neumannBoundary(), read_function_space=V, write_object=V,
fixed_boundary=clampedBoundary())
precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary)

fenics_dt = precice_dt # if fenics_dt == precice_dt, no subcycling is applied
# fenics_dt = 0.02 # if fenics_dt < precice_dt, subcycling is applied
dt = Constant(np.min([precice_dt, fenics_dt]))

# clamp the beam at the bottom
bc = DirichletBC(V, Constant((0, 0)), clampedBoundary())
bc = DirichletBC(V, Constant((0, 0)), fixed_boundary)

# alpha method parameters
alpha_m = Constant(0.2)
Expand Down
11 changes: 4 additions & 7 deletions FSI/flap_perp_2D/SU2-FEniCS/Solid/perp-flap.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,6 @@
from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \
TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, project, \
Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system
import dolfin

from ufl import nabla_div
import numpy as np
import matplotlib.pyplot as plt
Expand Down Expand Up @@ -63,26 +61,25 @@ def neumann_boundary(x, on_boundary):

# Function to calculate displacement Deltas
u_delta = Function(V)
u_ref = Function(V)

f_N_function = interpolate(Expression(("1", "0"), degree=1), V)
u_function = interpolate(Expression(("0", "0"), degree=1), V)

coupling_boundary = AutoSubDomain(neumann_boundary)
fixed_boundary = AutoSubDomain(clamped_boundary)

precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json")

clamped_boundary_domain = AutoSubDomain(clamped_boundary)
force_boundary = AutoSubDomain(neumann_boundary)

# Initialize the coupling interface
precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=clamped_boundary_domain)
precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary)

fenics_dt = precice_dt # if fenics_dt == precice_dt, no subcycling is applied
# fenics_dt = 0.02 # if fenics_dt < precice_dt, subcycling is applied
dt = Constant(np.min([precice_dt, fenics_dt]))

# clamp the beam at the bottom
bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)
bc = DirichletBC(V, Constant((0, 0)), fixed_boundary)

# alpha method parameters
alpha_m = Constant(0.2)
Expand Down