How to find border markings

1

Hello everyone.
I just started to study fenics.
I have a question. I want to set the boundary conditions as a class if the boundaries == 2 and boundaries == 3.
I tried to do it but it did not work.
Thanks in advance!

from dolfin import *
import matplotlib.pyplot as plt

mesh = Mesh("cube.xml")
mesh.init()
subdomains = MeshFunction("size_t", mesh, "cube_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "cube_facet_region.xml")

dx = Measure('dx', domain=mesh, subdomain_data = subdomains)
ds = Measure('ds', domain=mesh, subdomain_data = boundaries)

order = 1
BDM = FiniteElement("BDM", triangle, order)
DG = FiniteElement("DG", triangle, order - 1)
mixed = MixedElement(BDM, DG, DG)

W = FunctionSpace(mesh, mixed)

tdim = mesh.topology().dim()
mesh.init(tdim-1, tdim)

dw = TrialFunction(W)

wt = TestFunction(W)

v, q, r = split(wt)

w  = Function(W)
w0 = Function(W)

u , p, s = split(w)
u0, p0, s0 = split(w0)

sm = 0.5 * (s0 + s)

class PressureBC(UserExpression):
	def eval_cell(self, values, x, cell):
		for facet in facets(mesh):
			if facet.index() == 2:
				values[0] = 1

class SaturationBC(UserExpression):
	def eval(self, values, x):
		for facet in facets(mesh):
			if facet.index() == 3:
				values[0] = 1

pbar = PressureBC()
sbar = SaturationBC()


dt = Constant(0.01)

kinv = Expression("1 / max (exp(-pow(10 * (x[1] - 0.5 - 0.1 * sin(10 * x[0])), 2)), 0.01)", degree=2)
zero = Constant(0.0)
Kinv = as_matrix(((kinv, zero), (zero, kinv)))

mu_rel = Constant(0.2)

def linv(s):
	return 1.0 / ((1.0 / mu_rel)*s**2 + (1.0 - s)**2)

def F(s):
	return s**2 / (s**2 + mu_rel * (1.0 - s)**2)

n = FacetNormal(mesh)

un = 0.5 *( inner(u0, n) + sqrt(inner(u0, n) * inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - sqrt(inner(u0, n) * inner(u0, n)))
stab = inner(jump(r), un('+')*F(sm)('+')-un('-')*F(sm)('-'))*dS+r*un_h*sbar*ds

#stab = inner(jump(r), un('+')*F(sm)('+')-un('-')*F(sm)('-'))*ds+r*(un_h)*sbar*ds

F1 = inner(v, linv(sm)*Kinv*u)*dx - div(v)*p*dx + inner(v, pbar*n)*ds
F2 = q * div(u) * dx
F3 = r * (s - s0) * dx - dt * inner(grad(r), F(sm) * u) * dx + dt * r * F(sm) * un * ds + dt * stab

F = F1 + F2 + F3
J = derivative(F, w, dw)
info("Hello")
print(" ")

newton_solver_parameters = {"nonlinear_solver": "newton",
                          "newton_solver": {"linear_solver": "bicgstab", "preconditioner": "default", "absolute_tolerance" : 1e-6, "relative_tolerance" :1e-5, "maximum_iterations": 10}}

problem = NonlinearVariationalProblem(F, w, J=J)
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(newton_solver_parameters)
#info(solver.parameters, True)


# Solve the problem
(iter, converged) = solver.solve()

# Check for convergence
if not converged:
    warning("This demo is a complex nonlinear problem. Convergence is not guaranteed when modifying some parameters or using PETSC 3.2.")



folder = "res/examp/"

u_file = File(folder + "velocity.pvd")
p_file = File(folder + "pressure.pvd")
s_file = File(folder + "saturation.pvd")

max_time_iters = 100
for iter in range(max_time_iters + 1):
	print("Time iteration = %d" % iter)
	
	
	w0.assign(w)
	solver.solve()

	u, p, s = w.split()

	u_file << u
	p_file << p
	s_file << s
2

A couple of pointers. Please supply a minimal working example next time you ask a question. This should be something that is reproducable by anyone on any computer with dolfin installed.
As you are using a custom mesh, I cannot use the code you supplied.
Also, a minimal working example should be stripped of all code that is not required to get to the desired issue.

In your case, I would rather restrict the integration measure, instead of trying to restrict the expression (minimal example attached):

from dolfin import *

mesh = UnitCubeMesh(10,10,10)
mesh.init()
boundaries = MeshFunction("size_t", mesh, 2, 0)


class BoundaryMarker(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < 0.2 and on_boundary
BoundaryMarker().mark(boundaries, 1)

ds = Measure('ds', domain=mesh, subdomain_data = boundaries)

class PressureBC(UserExpression):
    def eval_cell(self, values, x, cell):
        values[0] = 1

pbar = PressureBC()
print("Restricted ", assemble(pbar*ds(1)))
print("Unrestricted ", assemble(pbar*ds))

Which gives you the output:

Restricted  1.3999999999999964
Unrestricted  5.9999999999999005
4

First thing i would do is to move this until after you have solved the problem.

5

I’m sorry I did not understand

6

You are assigning w0 to w before your are solving the PDE in your loop. This update should be done after you have solved the pde.

7

I get the same results. Did I spell classes correctly?

8

As you have not provided a minimal example illustrating your issue, it is a rather tedious thing to debug your code. As you have not supplied any reproducible meshes, or changed your example to illustrate how one should actually restrict the boundary conditions, there is limits to what I am willing to do by visual inspection of your code.

You should use the minimal example I showed above as a template to illustrate the issues you are currently having. For more guidelines see Read before posting: How do I get my question answered?

9

Thank you very much for your attention and help me.
In general, I have such a grid. And I want to set the pressure at border 2 and the saturation at border 3.

LX = 2571 * 2;
LY = 2571;
r = 0.1;
x1 = LX/2 - LX/4;
y1 = LY/2;

x2 = LX/2 + LX/4;
y2 = LY/2;

p1 = 500;
p2 = 300;

Point(1) = {0, 0, 0, p1};
Point(2) = {LX, 0, 0, p1};
Point(3) = {LX, LY, 0, p1};
Point(4) = {0, LY, 0, p1};

Point(5) = {x1, y1, 0, p2};
Point(6) = {x1 + r, y1, 0, p2};
Point(7) = {x1 - r, y1, 0, p2};
Point(8) = {x1, y1 + r, 0, p2};
Point(9) = {x1, y1 - r, 0, p2};

Point(10) = {x2, y2, 0, p2};
Point(11) = {x2 + r, y2, 0, p2};
Point(12) = {x2 - r, y2, 0, p2};
Point(13) = {x2, y2 + r, 0, p2};
Point(14) = {x2, y2 - r, 0, p2};

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};

Circle(5) = {6, 5, 9};
Circle(6) = {9, 5, 7};
Circle(7) = {7, 5, 8};
Circle(8) = {8, 5, 6};

Circle(9) = {11, 10, 14};
Circle(10) = {14, 10, 12};
Circle(11) = {12, 10, 13};
Circle(12) = {13, 10, 11};

Line Loop(13) = {4, 1, 2, 3};
Line Loop(14) = {12, 9, 10, 11};
Line Loop(15) = {7, 8, 5, 6};
Plane Surface(16) = {13, 14, 15};

Physical Line(1) = {3, 4, 1, 2};
Physical Line(2) = {8, 7, 6, 5, 11, 10, 9, 12};
Physical Line(3) = {4};
Physical Surface(1) = {16};

Plane Surface(17) = {15};
Plane Surface(18) = {14};
Physical Surface(2) = {18, 17};