Skip to content

Add periodic bcs #35

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Jan 22, 2021
Merged

Add periodic bcs #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
176fc0d
Extended GmshDiscreteModel to periodic bcs.
fverdugo Jan 21, 2021
24bca15
Merge branch 'master' of github.com:gridap/GridapGmsh.jl into add_per…
fverdugo Jan 21, 2021
bcb7deb
Periodic BCs working (tested in 2D)
fverdugo Jan 22, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 70 additions & 7 deletions src/GmshDiscreteModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,9 @@ function GmshDiscreteModel(mshfile; renumber=true)
renumber && gmsh.model.mesh.renumberElements()

grid, cell_to_entity = _setup_grid(gmsh)
grid_topology = UnstructuredGridTopology(grid)
labeling = _setup_labeling(gmsh,grid,grid_topology,cell_to_entity)
cell_to_vertices, vertex_to_node, node_to_vertex = _setup_vertices(gmsh,grid)
grid_topology = UnstructuredGridTopology(grid,cell_to_vertices,vertex_to_node)
labeling = _setup_labeling(gmsh,grid,grid_topology,cell_to_entity,vertex_to_node,node_to_vertex)

gmsh.finalize()

Expand Down Expand Up @@ -51,6 +52,53 @@ function _setup_grid(gmsh)

end

function _setup_vertices(gmsh,grid)
cell_to_nodes = get_cell_node_ids(grid)
dimTags = gmsh.model.getEntities()
if _has_periodic_bcs(gmsh,dimTags)
cell_to_vertices, vertex_to_node, node_to_vertex = _setup_vertices_periodic(
gmsh,dimTags,cell_to_nodes,num_nodes(grid))
else
cell_to_vertices = cell_to_nodes
vertex_to_node = 1:num_nodes(grid)
node_to_vertex = vertex_to_node
end
cell_to_vertices, vertex_to_node, node_to_vertex
end

function _has_periodic_bcs(gmsh,dimTags)
for (dim,tag) in dimTags
tagMaster, nodeTags, = gmsh.model.mesh.getPeriodicNodes(dim,tag)
if length(nodeTags) > 0
return true
end
end
return false
end

function _setup_vertices_periodic(gmsh,dimTags,cell_to_nodes,nnodes)
# Assumes linear grid
Copy link

@stevengj stevengj Jan 29, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What does this mean?

Copy link
Member Author

@fverdugo fverdugo Jan 29, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

High order meshes not yet supported, but it should be doable as well. Independently of this, you can use high order fe spaces in gridap as always.

node_to_node_master = fill(UNSET,nnodes)
_node_to_node_master!(node_to_node_master,gmsh,dimTags)
slave_to_node_slave = findall(node_to_node_master .!= UNSET)
slave_to_node_master = node_to_node_master[slave_to_node_slave]
node_to_vertex = fill(UNSET,nnodes)
vertex_to_node = findall(node_to_node_master .== UNSET)
node_to_vertex[vertex_to_node] = 1:length(vertex_to_node)
node_to_vertex[slave_to_node_slave] = node_to_vertex[slave_to_node_master]
cell_to_vertices = Table(lazy_map(Broadcasting(Reindex(node_to_vertex)),cell_to_nodes))
cell_to_vertices, vertex_to_node, node_to_vertex
end

function _node_to_node_master!(node_to_node_master,gmsh,dimTags)
for (dim,tag) in dimTags
tagMaster, nodeTags, nodeTagsMaster, = gmsh.model.mesh.getPeriodicNodes(dim,tag)
if length(nodeTags) > 0
node_to_node_master[nodeTags] .= nodeTagsMaster
end
end
end

function _setup_dim(gmsh)
entities = gmsh.model.getEntities()
D = -1
Expand Down Expand Up @@ -298,7 +346,7 @@ function _fill_cell_to_entity!(cell_to_entity,elemTags,entity,o)
end
end

function _setup_labeling(gmsh,grid,grid_topology,cell_to_entity)
function _setup_labeling(gmsh,grid,grid_topology,cell_to_entity,vertex_to_node,node_to_vertex)

D = num_cell_dims(grid)
dim_to_gface_to_nodes, dim_gface_to_entity = _setup_faces(gmsh,D)
Expand All @@ -313,6 +361,8 @@ function _setup_labeling(gmsh,grid,grid_topology,cell_to_entity)
node_to_label = _setup_node_to_label(gmsh,dim_to_offset,nnodes)

labeling = _setup_face_labels(
vertex_to_node,
node_to_vertex,
grid_topology,
dim_to_gface_to_nodes,
dim_gface_to_entity,
Expand Down Expand Up @@ -405,6 +455,8 @@ function _setup_node_to_label(gmsh,dim_to_offset,nnodes)
end

function _setup_face_labels(
vertex_to_node,
node_to_vertex,
grid_topology,
dim_to_gface_to_nodes,
dim_gface_to_entity,
Expand All @@ -413,16 +465,19 @@ function _setup_face_labels(
dim_to_group_to_name,
node_to_label)

vertex_to_label = node_to_label[vertex_to_node]

D = num_cell_dims(grid_topology)

dim_to_face_to_label = [node_to_label,]
dim_to_face_to_label = [vertex_to_label,]

for d in 1:D
z = fill(UNSET,num_faces(grid_topology,d))
push!(dim_to_face_to_label,z)
end

_fill_dim_to_face_to_label!(
node_to_vertex,
dim_to_face_to_label,
grid_topology,
dim_to_gface_to_nodes,
Expand All @@ -439,6 +494,7 @@ function _setup_face_labels(
end

function _fill_dim_to_face_to_label!(
node_to_vertex,
dim_to_face_to_label,
grid_topology,
dim_to_gface_to_nodes,
Expand All @@ -458,6 +514,7 @@ function _fill_dim_to_face_to_label!(
face_to_nodes = get_faces(grid_topology,d,0)
node_to_faces = get_faces(grid_topology,0,d)
gface_to_face = _setup_gface_to_face(
node_to_vertex,
face_to_nodes,
node_to_faces,
gface_to_nodes)
Expand Down Expand Up @@ -497,11 +554,13 @@ function _apply_offset_for_faces!(face_to_label,gface_to_entity,gface_to_face,of
end

function _setup_gface_to_face(
node_to_vertex,
face_to_nodes,
node_to_faces,
gface_to_nodes)

gface_to_face = _find_gface_to_face(
node_to_vertex,
face_to_nodes.data,
face_to_nodes.ptrs,
node_to_faces.data,
Expand All @@ -514,6 +573,7 @@ function _setup_gface_to_face(
end

function _find_gface_to_face(
node_to_vertex,
face_to_nodes_data,
face_to_nodes_ptrs,
node_to_faces_data::AbstractVector{T},
Expand All @@ -522,13 +582,14 @@ function _find_gface_to_face(
gface_to_nodes_ptrs) where T

ngfaces = length(gface_to_nodes_ptrs) - 1
gface_to_face = zeros(T,ngfaces)
gface_to_face = fill(T(UNSET),ngfaces)
n = max_cells_arround_vertex(node_to_faces_ptrs)
faces_around = fill(UNSET,n)
faces_around_scratch = fill(UNSET,n)

_fill_gface_to_face!(
gface_to_face,
node_to_vertex,
face_to_nodes_data,
face_to_nodes_ptrs,
node_to_faces_data,
Expand All @@ -544,6 +605,7 @@ end

function _fill_gface_to_face!(
gface_to_face,
node_to_vertex,
face_to_nodes_data,
face_to_nodes_ptrs,
node_to_faces_data,
Expand All @@ -566,16 +628,17 @@ function _fill_gface_to_face!(

for lnode in 1:nlnodes
node = gface_to_nodes_data[lnode+a]
vertex = node_to_vertex[node]
if lnode == 1
nfaces_around = _fill_cells_around_scratch!(
faces_around,
node,
vertex,
node_to_faces_data,
node_to_faces_ptrs)
else
nfaces_around_scratch = _fill_cells_around_scratch!(
faces_around_scratch,
node,
vertex,
node_to_faces_data,
node_to_faces_ptrs)
_set_intersection!(
Expand Down
19 changes: 19 additions & 0 deletions test/GmshDiscreteModelsTests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,4 +50,23 @@ model = GmshDiscreteModel(mshfile)
test_discrete_model(model)
check_interpolation(model)

mshfile = joinpath(@__DIR__,"periodic.msh")
model = GmshDiscreteModel(mshfile)
test_discrete_model(model)
u(x) = 2*x[1]
reffe = ReferenceFE(lagrangian,Float64,1)
V = TestFESpace(model,reffe,dirichlet_tags=["L","PL"])
U = TrialFESpace(V,u)
Ω = Triangulation(model)
Γ = BoundaryTriangulation(model,tags="R")
n_Γ = get_normal_vector(Γ)
dΩ = Measure(Ω,2)
dΓ = Measure(Γ,2)
a(u,v) = ∫( ∇(u)⋅∇(v) )dΩ
l(v) = ∫( n_Γ⋅(v*∇(u)) )dΓ
op = AffineFEOperator(a,l,U,V)
uh = solve(op)
#writevtk(Ω,"sol",cellfields=["uh"=>uh,"u"=>u,"e"=>u-uh,"zh"=>zero(U)])
@test sum( ∫(abs2(u - uh))dΩ ) < 1.0e-9

end # module
40 changes: 40 additions & 0 deletions test/periodic.geo
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
//+
Point(1) = {-1, -1, 0, 0.25};
//+
Point(2) = {0, -1, 0, 2.0};
//+
Point(3) = {0, -0, 0, 0.25};
//+
Point(4) = {-1, 0, 0, 2.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 1};
//+
Curve Loop(1) = {4, 1, 2, 3};
//+
Plane Surface(1) = {1};

Periodic Line {3} = {1} Translate {0, 1, 0};
//+
Physical Point("PL") = {1};
//+
Physical Point("PR") = {2};
//+
Physical Point("PR") += {3};
//+
Physical Point("PL") += {4};
//+
Physical Curve("P") = {1};
//+
Physical Curve("R") = {2};
//+
Physical Curve("P") += {3};
//+
Physical Curve("L") = {4};
//+
Physical Surface("S") = {1};
Loading