I’ve changed the last segment by deliting if facet.exterior() statement and the snipped worked fine!
boundaries.set_all(0)
cell_to_facet = mesh.topology()(mesh.topology().dim(), mesh.topology().dim()-1)
for cell in cells(mesh):
for face in cell_to_facet(cell.index()):
facet = Facet(mesh, face)
if facet.exterior():
if domains.array()[cell.index()]==1:
boundaries.set_value(facet.index(), 1)
It looks like it worked. Thank you for your help!
Hello @dokken. I have another question that is similar to the previous one.
I have successfully created boundary condition to use in solving poisson equation.
Now I want to improve the code so it would slove poisson equation with time stepping.
I chose the formulation from the example:
F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
Now I have to create the initial state - the u_n.
In the example the u_n is creating via interpolation:
u_D = Constant(0)
u_n = interpolate(u_D, V)
How can I convert meshFunction boundaries defined with your snipped from this answer
I’ve tried this:
u_n = project(boundaries, V)
But obviously failed with the following:
Traceback (most recent call last):
File "fenics/poisson_time.py", line 51, in <module>
u_n = project(boundaries, V)
File "/usr/lib/python3/dist-packages/dolfin/fem/projection.py", line 129, in project
L = ufl.inner(w, v) * dx
File "/usr/lib/python3/dist-packages/ufl/operators.py", line 166, in inner
b = as_ufl(b)
File "/usr/lib/python3/dist-packages/ufl/constantvalue.py", line 472, in as_ufl
" to any UFL type." % str(expression))
ufl.log.UFLValueError: Invalid type conversion: <dolfin.cpp.mesh.MeshFunctionSizet object at 0x7f33265fbd70> can not be converted to any UFL type.
Hi @mkondratyev85
First of all, u_n has to be a Discontinuous Galerkin function of order 0. Then you can do the following:
from dolfin import *
mesh = UnitCubeMesh(10,10,10)
mesh.init()
class top_half(SubDomain):
def inside(self, x, on_boundary):
return (x[0]<=0.5)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
top_half().mark(domains,1)
V = FunctionSpace(mesh, "DG", 0)
v = Function(V)
v.vector()[:] = domains.array()
XDMFFile("bdfunc.xdmf").write(v)
Thank you for the answer. Unfortunately, I don’t understand what is the equivalent of u_n in the snipped provided. I have my domains defined like this:
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("/tmp/fenics/mf.xdmf") as infile:
infile.read(mvc, "bc")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
domains = mf
Then, after I try this:
V = FunctionSpace(mesh, "DG", 0)
u_n = Function(V)
u_n.vector()[:] = domains.array()
I get the following error:
Traceback (most recent call last):
File "fenics/poisson_time.py", line 75, in <module>
u_n.vector()[:] = domains.array()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to set local values of PETSc vector.
*** Reason: Size of values array is not equal to local vector size.
*** Where: This error was encountered inside PETScVector.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
My new question is: Is “bc” a field on a tetrahedron or a triangle (assuming the mesh consists of tetrahedron cells). If bc is on a tetrahedron, you should have a mesh value collection of the the same dimesion as the topology, and not dim-1. Then the array of values in mf should be of the same length as u_n.vector. This can be confirmed by printing len(mf.array()) and len(u_n.vector().array()).
I get the following error:
Traceback (most recent call last):
File "fenics/poisson_time.py", line 76, in <module>
print(len(u_n.vector().array()))
AttributeError: 'dolfin.cpp.la.PETScVector' object has no attribute 'array'
sorry, u_n.vector().get_local()
mf.array() -> 422400
u_n.vector().get_local() - > 78674
So what kind of function is the function that you are trying to load into dolfin, is it defined per cell or per facet?
What do you mean by function?
I try to load two meshes like these:
# Load mesh to solve
mesh = Mesh()
with XDMFFile("/tmp/fenics/mesh.xdmf") as infile:
infile.read(mesh)
# Load mesh with boundary conditions
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("/tmp/fenics/mf.xdmf") as infile:
infile.read(mvc, "bc")
domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)
And then I create another function like this:
# Convert cell values to face values in boundary condition mesh
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
cell_to_facet = mesh.topology()(mesh.topology().dim(), mesh.topology().dim()-1)
for cell in cells(mesh):
for face in cell_to_facet(cell.index()):
facet = Facet(mesh, face)
# if facet.exterior():
if domains.array()[cell.index()]==2:
boundaries.set_value(facet.index(), 2)
And here is the part where I try to define initial conditions based on loaded meshes:
V = FunctionSpace(mesh, 'P', 1)
# Define initial boundary conditions
u_n = Function(V)
print(len(mf.array()))
print(len(u_n.vector().get_local()))
u_n.vector()[:] = domains.array()
So you are trying to make discontinuous data that is defined per cell match up with a continuous functions space with degrees of freedom on each vertex? What should the value between to cells Which different internal values be?
It could be just an average of two cells.
It looks like this approach solves my problem:
u_D = Function(V)
dm = V.dofmap()
for cell in cells(mesh):
if domains[cell] < 2:
u_D.vector()[dm.cell_dofs(cell.index())] = 0
else:
u_D.vector()[dm.cell_dofs(cell.index())] = 100
u_n = interpolate(u_D, V)
But it looks very clumsy and unidiomatic. I’m sure there must be a better, more elegant way.
So if you really want it to be the data that you have in the file, the first function u_n has to ve a DG 0 function, as in my post above. Then after the first dolve you could update u_n to be a CG one function. If you really want it to be a CG 1 function at the first solve you can first load the data into a DG0 function, then project this to a CG1 function. I can Give an example of this next time Im at my computer
I would highly appreciate if you could.
Here you go:
from dolfin import *
mesh = UnitCubeMesh(10,10,10)
mesh.init()
class top_half(SubDomain):
def inside(self, x, on_boundary):
return (x[0]<=0.25)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
top_half().mark(domains,1)
V = FunctionSpace(mesh, "DG", 0)
v = Function(V)
v.vector()[:] = domains.array()
XDMFFile("bdfunc.xdmf").write(v)
W = FunctionSpace(mesh, "CG", 1)
u_n = interpolate(v, W)
print(u_n.vector().get_local())
XDMFFile("u_n.xdmf").write(u_n)
So as you can see, domains is a function per cell, Im moving it into a DG0 function, and then projecting it to the CG1 space.
W = FunctionSpace(mesh, "CG", 1)
Thank you so much for the help, @dokken! The method works!
I see you always save the output to the xdmf format. What advantages does it have compared to pvd format?
XDMF is the preferred format for parallel handling, and can be used to do checkpointing and can visualize other elements than CG 1, such as CG3 and DG elements.