Commit 9ce434cc authored by Tilman Aleman's avatar Tilman Aleman
Browse files

trying vectorlaplace in julia

parent 9cae2963
using Gridap
using GridapEmbedded
using LinearAlgebra
using GridapGmsh: gmsh
using GridapGmsh
# Define geometry
gmsh.initialize()
gmsh.clear()
box = gmsh.model.occ.addBox(
-2.,-2.,-2.,2.,2.,2.)
gmsh.model.occ.synchronize()
gmsh.model.mesh.setSize(gmsh.model.getEntities(0),0.1)
gmsh.model.mesh.generate(3)
gmsh.write("model.msh")
gmsh.finalize()
n = 10
# Background model
partition = (n,n,n)
pmin = 2*Point(-1,-1,-1)
pmax = 2*Point(1,1,1)
# bgmodel = simplexify(CartesianDiscreteModel(pmin,pmax,partition))
bgmodel = GmshDiscreteModel("model.msh")
# Select geometry
R = 0.5
ϕ(x) = norm(x)-1
# Forcing data
ud = 1
f = 10
println("doing cutter")
lscutter = GridapEmbedded.LevelSetCutters.sphere(1)
# Cut the background model
cutgeo = cut(bgmodel, lscutter)
# Setup integration meshes
Ω = Triangulation(cutgeo,PHYSICAL)
Γd = EmbeddedBoundary(cutgeo)
# Setup normal vectors
n_Γd = get_normal_vector(Γd)
#writevtk(Ω,"trian_O")
# writevtk(Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
#writevtk(Γg,"trian_Gg",cellfields=["normal"=>n_Γg])
#writevtk(Triangulation(bgmodel),"bgtrian")
# Setup Lebesgue measures
order = 2
degree = 2*order
= Measure(Ω,degree)
dΓd = Measure(Γd,degree)
reffe = ReferenceFE(lagrangian, VectorValue{3, Float64}, order)
# Setup FESpace
Ω_act = Triangulation(cutgeo,ACTIVE)
normalspace = FESpace(Ω_act,reffe)
# ntilde = interpolate_everywhere(n_Γd, normalspace)
println("setup TnT")
V = TestFESpace(Ω_act,ReferenceFE(lagrangian,VectorValue{3,Float64},order),conformity=:H1)
U = TrialFESpace(V)
# Weak form
nex = (ϕ)
h = (pmax - pmin)[1] / partition[1]
Pmats(n) = (x)->identity(x)-(n⊗ n) x
Pmats2(n) = (x)-> identity(x)-x⋅ (n⊗n)
Pmat = Pmats(n_Γd)
Pmat2 = Pmats2(n_Γd)
ϵ = 1
η = 1.0/h^2
ρu = 1.0/h
println("setup blf")
E(u) = Pmat2(Pmat(((u)+transpose((u)))))
a1(u,v)= ( E(u) E(v) )*dΓd
a2(u,v) = ( ϵ*((Pmat(u))(Pmat(v))))*dΓd
a3(u,v) = (η*((u n_Γd)(v⋅n_Γd)))*dΓd
# a4(u,v) = ∫(ρu*((∇(u)⋅n_Γd)⋅(∇(v)⋅n_Γd)))*dΩ
a4(u,v) = (ρu*(((u)⋅nex)((v)⋅nex)))*
a(u,v) =(a1(u,v)+a2(u,v)+a3(u,v)+a4(u,v))
println("setup lf")
l(v) =(VectorValue{3,Float64}(1,1,1) v)*dΓd
# ∫( v*f ) * dΩ +
# ∫( (γd/h)*v*ud - (n_Γd⋅∇(v))*ud ) * dΓd
# FE problem
println("setup FEop")
op = AffineFEOperator(a,l,U,V)
@time AffineFEOperator(a,l,U,V)
println("solve")
ls = LinearFESolver()
uh = solve(op)
# Postprocess
# if outputfile !== nothing
writevtk(Ω,"vectorlaplace_julia",cellfields=["uh"=>uh])
# end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment