diff --git a/FSI/flap_perp_2D/OpenFOAM-FEniCS/Solid/perp-flap.py b/FSI/flap_perp_2D/OpenFOAM-FEniCS/Solid/perp-flap.py index 4b1891d17..ee19aefca 100644 --- a/FSI/flap_perp_2D/OpenFOAM-FEniCS/Solid/perp-flap.py +++ b/FSI/flap_perp_2D/OpenFOAM-FEniCS/Solid/perp-flap.py @@ -1,9 +1,7 @@ # 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 @@ -11,22 +9,17 @@ 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 @@ -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) @@ -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) diff --git a/FSI/flap_perp_2D/SU2-FEniCS/Solid/perp-flap.py b/FSI/flap_perp_2D/SU2-FEniCS/Solid/perp-flap.py index 5376588a9..a57f922a3 100644 --- a/FSI/flap_perp_2D/SU2-FEniCS/Solid/perp-flap.py +++ b/FSI/flap_perp_2D/SU2-FEniCS/Solid/perp-flap.py @@ -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 @@ -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)