Try again now (i’ve updated my code once again), you simply need to replace the " with properly formatted once.
1 Like
Hi dokken! Thanks for the solution.
I have one more question:
Does “mf” preserve the physical tags I attached to the boundaries in gmsh?
I have this question because when I’m executing the code with n=FacetNormal(mesh), it is not recognizing the boundaries on which “n” should be defined.
This is dependent of how you have written your Gmsh file. Try for instance to open your “mf.xdmf” in Paraview and look at what boundaries has been saved.
Trying to open mf.xdmf in Paraview (version 5.9.0) using XDMF Reader is giving me the error message “Failed to read data”.
I tried that. But everytime I try it, Paraview crashes.
Here is my gmsh script for the structure without the mesh:
//+
SetFactory("OpenCASCADE");
Circle(1) = {0, 0, 0, 5, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 3, 0, 2*Pi};
//+
Line Loop(1) = {1};
//+
Line Loop(2) = {2};
//+
Plane Surface(1) = {1, 2};
//+
Physical Surface("S1", 3) = {1};
//+
Physical Line("C1", 1) = {1};
//+
Physical Line("C2", 2) = {2};Why are you saying that your subdomain id is 12? In the mesh
you never use the tag 12, but 1,2,3.
I have tested your mesh with gmsh v 4.6.0 and ran the script I supplied above with no issues, only with a single change:
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)
print(assemble(1*ds_custom(1)) ,assemble(1*ds_custom(2)))
This works perfectly fine and also preserves the structure from gmsh. But I think with my FEniCS script, the problem is with the code n=FacetNormal(mesh) . Somehow, it is not defining the normal to the boundaries properly. I need to figure that out.
What do you mean by it doesn’t define the normal properly?
Could you produce a minimal code example that illustrates what you are talking about?
Before I post the code, I need to make sure of one more thing.
If I want to impose a Dirichlet BC (u=0) on the outer circle (physical tag: 1), will the following code work?
c=Expression(("0.0", "0.0"),degree=2)
uout=DirichletBC(V, c, 1)
bc=[uout]No, you need to supply the mesh tag as well, i.e.:
c=Expression(("0.0", "0.0"),degree=2)
uout=DirichletBC(V, c, mf, 1)
bc=[uout]
1 Like
Thanks dokken! Problem solved. Without the essential BC, the solution to my problem was not a unique one. So it was showing a very small translation. As a result, the 3D Glyphs were all oriented the wrong way.
Your patience and sharing your knowledge and experience with someone who is new to FEniCS is really appreciated. Thanks a lot.
As this is a rather long question, with a lot of information in it that does not directly relate the the original post, I would prefer if you deleted this post and made it a new question. This way we avoid long posts where the solution to a problem becomes unclear due to the amount of posts.
Dear sir’
used above syntax to import gmsh file to fenics. but i am not able to define function on mesh . i am attaching code and error
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data(“gmsh:physical”, cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={“name_to_read”:[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
import meshio
msh = meshio.read(“subbygmsh.msh”)
triangle_mesh = create_mesh(msh, “triangle”, True)
line_mesh = create_mesh(msh, “line”, True)
meshio.write(“mesh1.xdmf”, triangle_mesh)
meshio.write(“mf1.xdmf”, line_mesh)
from dolfin import*
mesh=Mesh()
mvc = MeshValueCollection(“size_t”, mesh, mesh.topology().dim())
with XDMFFile(“mesh1.xdmf”) as infile:
infile.read(mesh)
infile.read(mvc, “name_to_read”)
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds= Measure(“dx”, domain=mesh, subdomain_data=cf)
mvc = MeshValueCollection(“size_t”, mesh, mesh.topology().dim()-1)
with XDMFFile(“mf1.xdmf”) as infile:
infile.read(mvc, “name_to_read”)
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds= Measure(“ds”, domain=mesh, subdomain_data=mf)
V=FunctionSpace(mesh, “P”, 1)
u=TrialFunction(V)
v=TestFunction(V)
f=Constant(1.0)
u0=Constant(0.0)
bc=DirichletBC(V, u0, mf, 3)
a0=1
a1=2
A=a0*inner(grad(u), grad(v))dx(1)+a1inner(grad(u), grad(v))*dx(2)
L=fvds(4)
u=Function(V)
and error is
*** The MPI_Comm_rank() function was called after MPI_FINALIZE was invoked.
*** This is disallowed by the MPI standard.
*** Your MPI job will now abort.
[gangadhar-Inspiron-3593:24187] Local abort after MPI_FINALIZE started completed successfully, but am not able to aggregate error messages, and not able to guarantee that all other processes were killed!
Please format your code using ``` and make sure it is possible to copy-paste the code into an empty editor, save it and run it.
Currently, your example is not reproducible as you have not supplied your domain.
Finally, please report which line your error occurs in.
thank you for your replays . and i am learning lots of thing from your replays.
this is .geo file for domain.
// Gmsh project created on Fri Mar 5 14:40:36 2021
SetFactory(“OpenCASCADE”);
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {1, 0, 0, 1.0};
//+
Point(3) = {0, 1, 0, 1.0};
//+
Point(4) = {1, 1, 0, 1.0};
//+
Point(5) = {0, 0.5, 0, 1.0};
//+
Point(6) = {1, 0.5, 0, 1.0};
//+
Line(1) = {3, 5};
//+
Line(2) = {3, 4};
//+
Line(3) = {4, 6};
//+
Line(4) = {5, 6};
//+
Line(5) = {5, 1};
//+
Line(6) = {6, 2};
//+
Line(7) = {2, 1};
//+
Curve Loop(1) = {1, 4, -3, -2};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {7, -5, 4, 6};
//+
Plane Surface(2) = {2};
//+
Physical Surface(1) = {1};
//+
Physical Surface(2) = {2};
//+
Physical Curve(3) = {2};
//+
Physical Curve(4) = {3, 6};
and this is fenics code
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data(“gmsh:physical”, cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={“name_to_read”:[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
import meshio
msh = meshio.read(“subbygmsh.msh”)
triangle_mesh = create_mesh(msh, “triangle”, True)
line_mesh = create_mesh(msh, “line”, True)
meshio.write(“mesh1.xdmf”, triangle_mesh)
meshio.write(“mf1.xdmf”, line_mesh)
from dolfin import*
mesh=Mesh()
mvc = MeshValueCollection(“size_t”, mesh, mesh.topology().dim())
with XDMFFile(“mesh1.xdmf”) as infile:
infile.read(mesh)
infile.read(mvc, “name_to_read”)
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds= Measure(“dx”, domain=mesh, subdomain_data=cf)
mvc = MeshValueCollection(“size_t”, mesh, mesh.topology().dim()-1)
with XDMFFile(“mf1.xdmf”) as infile:
infile.read(mvc, “name_to_read”)
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds= Measure(“ds”, domain=mesh, subdomain_data=mf)
V=FunctionSpace(mesh, “P”, 1)
u=TrialFunction(V)
v=TestFunction(V)
f=Constant(1.0)
u0=Constant(0.0)
bc=DirichletBC(V, u0, mf, 3)
a0=1
a1=2
A=a0*inner(grad(u), grad(v))dx(1)+a1inner(grad(u), grad(v))*dx(2)
L=fvds(4)
u=Function(V)
and error is generated by u=Function(V)
and error is
*** The MPI_Comm_rank() function was called after MPI_FINALIZE was invoked.
*** This is disallowed by the MPI standard.
*** Your MPI job will now abort.
[gangadhar-Inspiron-3593:24187] Local abort after MPI_FINALIZE started completed successfully, but am not able to aggregate error messages, and not able to guarantee that all other processes were killed!
Please format your code using ``` encapsulation, as I mentioned in the previous post. Otherwise it is not possible to copy paste the code without tons of formatting
sorry sir , am not getting what you’re telling but i tried . please help me
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
import meshio
msh = meshio.read("subbygmsh.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh1.xdmf", triangle_mesh)
meshio.write("mf1.xdmf", line_mesh)
from dolfin import*
mesh=Mesh()
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mesh1.xdmf") as infile:
infile.read(mesh)
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds= Measure("dx", domain=mesh, subdomain_data=cf)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf1.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds= Measure("ds", domain=mesh, subdomain_data=mf)
V=FunctionSpace(mesh, "P", 1)
u=TrialFunction(V)
v=TestFunction(V)
f=Constant(1.0)
u0=Constant(0.0)
bc=DirichletBC(V, u0, mf, 3)
a0=1
a1=2
A=a0*inner(grad(u), grad(v))*dx(1)+a1*inner(grad(u), grad(v))*dx(2)
L=f*v*ds(4)
u=Function(V)
I cannot reproduce your issue in a docker container:
docker run -ti -v $PWD:/home/shared -w /home/shared --rm quay.io/fenicsproject/dev
export HDF5_DIR="/usr/lib/x86_64-linux-gnu/hdf5/mpich"
export HDF5_MPI="ON"
export CC=mpicc
pip3 install gmsh --user
pip3 install --no-binary=h5py h5py meshio --user
sudo apt-get update
sudo apt-get install -y -qq python3-pip libglu1 libxrender1 libxft2 libxinerama1
sudo apt-get install libxcursor1
/home/fenics/.local/lib/python3.6/site-packages/gmsh-4.8.0-Linux64-sdk/bin/gmsh -2 subbygmsh.geo
How did you install dolfin and what kind of system are you running on?
i am using ubuntu and i installed just fenics by sudo apt-get install fenics