Hello,
I have a problem with applying multiple constant b.c. in Fenicsx (fenics-dolfinx 0.7.0.dev0). I have inspired from this example: Setting multiple Dirichlet, Neumann, and Robin conditions — FEniCSx tutorial
I am trying to simulate a 2D plate that is being subjected to tension. The bottom and top ends have Neumman boundary conditions, and I have also added a node in the top right corner of the plate with a Dirichlet boundary condition to ensure stability but the uh is obtained as ([inf, inf, inf, inf, …]).
I think, I made an error in assigning the Dirichlet boundary condition to the node. Could you please let me know if you have any suggestions?
I used to work with fenics and mshr and I was able to apply such BC as:
class BoundaryC1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1, 0.01) and near(x[1], 1, 0.01)
Thanks!
Here is the code for the example that I am attempting to solve:
import numpy as np
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, div, dot, dx, grad, inner, lhs, rhs)
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
x = SpatialCoordinate(mesh)
# Define physical parameters
f=Constant(mesh, ScalarType(0))
kappa = Constant(mesh, ScalarType(1))
# Define function space and standard part of variational form
V = FunctionSpace(mesh, ("CG", 1))
u, v = TrialFunction(V), TestFunction(V)
F = kappa * inner(grad(u), grad(v)) * dx - inner(f, v) * dx
boundaries = [(1, lambda x: np.isclose(x[1], 0)),
(2, lambda x: np.isclose(x[1], 1)),
(3, lambda x: np.isclose(x[0],1) & np.isclose(x[1],1).all())]
g1=Constant(mesh, ScalarType(-1))
g2=Constant(mesh, ScalarType(-1))
u1= lambda x: x[1]+x[0]-2
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = facet_tag.find(marker)
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
b1 = Constant(mesh, ScalarType(0))
# Define the BC
boundary_conditions = [BoundaryCondition("Neumann", 1, g1),
BoundaryCondition("Neumann", 2, g2),
BoundaryCondition("Dirichlet", 3, u1)]
bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
