Here is the .geo file of my mesh.
// Gmsh project created on Wed Jun 16 11:22:49 2021
SetFactory("OpenCASCADE");
//+
Point(1) = {0, 0, 0, 3.0};
//+
Point(2) = {2, 0, 0, 3.0};
//+
Point(3) = {2, 1, 0, 3.0};
//+
Point(4) = {0, 1, 0, 3.0};
//+
Line(1) = {4, 1};
//+
Line(2) = {1, 2};
//+
Line(3) = {2, 3};
//+
Line(4) = {3, 4};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Physical Curve("fixed", 5) = {1};
//+
Physical Curve("free", 6) = {3};
//+
Physical Curve("top", 7) = {4};
//+
Physical Curve("bottom", 8) = {2};
//+
Physical Surface("fluid", 9) = {1};
This is the code which I have written for built-in Unitsquare mesh.
from fenics import *
import matplotlib.pyplot as plt
# Problem parameters
L = 1
F = 10
EA = 2
p = F/EA
# Defining mesh
mesh = UnitSquareMesh(2, 1)
V = FunctionSpace(mesh, 'P', 1)
plot(mesh)
# Define boundary segments for Neumann, Robin and Dirichlet conditions
# Create mesh function over cell facets
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
# Mark lower boundary facets as subdomain 0
class LowerNeumanBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[1]) < tol
Gamma_R = LowerNeumanBoundary()
Gamma_R.mark(boundary_parts, 0)
# Mark upper boundary facets as subdomain 1
class UpperNeumannBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[1] - 1) < tol
Gamma_N = UpperNeumannBoundary()
Gamma_N.mark(boundary_parts, 1)
# Mark left boundary as subdomain 2
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0]) < tol
Gamma_0 = LeftBoundary()
Gamma_0.mark(boundary_parts, 2)
# Mark right boundary as subdomain 3
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0] - 1) < tol
Gamma_1 = RightBoundary()
Gamma_1.mark(boundary_parts, 3)
bc = DirichletBC(V, Constant(0), boundary_parts, 2)
ds = Measure("ds", subdomain_data=boundary_parts)
# Defining the problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = -inner(grad(u), grad(v))*dx
g = Expression('p', p=p, degree=1)
g1 = Constant(0)
L = f*v*dx - g*v*ds(3) + g1*v*ds(0) + g1*v*ds(1)
#problem = VariationalProblem(a, L, bc)
#u = problem.solve()
u = Function(V)
solve(a == L, u, bc)
plot(u)
u.vector().get_local()
array([ 0. , 0. , 2.5, 2.5, 5. , 5. ])
This is the code that I want to write for the mesh created in GMSH, and want to get the similar result of what got in the above code.
from fenics import *
import matplotlib.pyplot as plt
pip install --user meshio
import meshio
mesh_from_file = meshio.read("loaded_bar.msh")
pip install --user h5py
import h5py
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
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
plot(mesh)
# Problem parameters
L = 2
F = 10
EA = 2
p = F/EA
# Define boundary segments for Neumann, Robin and Dirichlet conditions
# Create mesh function over cell facets
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
AttributeError Traceback (most recent call last)
<ipython-input-13-2980ce74ea92> in <module>
2
3 # Create mesh function over cell facets
----> 4 boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
AttributeError: module 'dolfin.mesh' has no attribute 'topology'
I got the above error.
