JordiManyer on refined-discrete-models
Added get_refined_polytopes (compare)
JordiManyer on refined-discrete-models
Compressed is_change_possible a… (compare)
JordiManyer on refined-discrete-models
Small change to a docstring (compare)
JordiManyer on refined-discrete-models
Compacted is_related code (compare)
JordiManyer on refined-discrete-models
Added check in AdaptedTriangula… (compare)
JordiManyer on refined-discrete-models
Added a new adapt method for fu… (compare)
norm(T) < some_number ? some_value : other_value
, however I get the error: no method matching isless(::Gridap.CellData.OperationCellField{ReferenceDomain}, ::Float64)
. Is there any other way to implement something like this? Any tutorial or test that I can look at for something similar::DistributedCellfield
at any point on the mesh. For example, println(evaluate(uh.fields.part,get_cell_coordinates(Ω.trians.part)[1][2]))
, returns an errorError:
│ exception =
│ AssertionError: A check failed
CartesianDiscreteModel()
this works fine, but when I try to evaluate on a un-structuredDiscretemodel()
it shows A check failed
error.
I'm still trying to understand how to set up the transient elastic response of a simple beam. I've tried quite a few different things for formulating the TransientAffineFEOperator
but I'm still getting errors. I've played around with it for a few days now and get different errors for different things I try. I'm quite confident that this is a problem with me and not Gridap.
Here is my code. Any help would be appreciated!
using Gridap
using Gridap.Geometry
const E = 70.0e9
const ν = 0.33
const λ = (E * ν) / ((1 + ν) * (1 - 2 * ν))
const μ = E / (2 * (1 + ν))
const density = 2710 # kg/m^3
σ(ε) = λ * tr(ε) * one(ε) + 2 * μ * ε
# for making the FEA grid.
L = 7
W = 2
H = 3
mult = 3
# ODE Solver config
Δt = 0.01
θ = 0.5
time_start = 0.0
time_end = 3.0
function make_model()
model = CartesianDiscreteModel((0, W, 0, L, 0, H), (W * mult, L * mult, H * mult))
labels = get_face_labeling(model)
add_tag_from_tags!(labels, "fixed", ["tag_24",
"tag_03",
"tag_04",
"tag_07",
"tag_08",
"tag_10",
"tag_12",
"tag_19",
"tag_20"
])
add_tag_from_tags!(labels, "load", ["tag_23"])
writevtk(model, "model_with_fixed")
model
end
function make_U0_solution(model; order=1)
reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
V0 = TestFESpace(model, reffe;
conformity=:H1,
dirichlet_tags=["fixed", "load"],
dirichlet_masks=[(true, true, true), (true, false, false)])
g2_disp(x) = VectorValue(0.05, 0.0, 0.0)
g1_fixed(x) = VectorValue(0.0, 0.0, 0.0)
# From functions `g1` and `g2`, we define the trial space as follows:
U = TrialFESpace(V0, [g1_fixed, g2_disp])
degree = 2 * order
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)
a(u, v) = ∫(ε(v) ⊙ (σ ∘ ε(u))) * dΩ
l(v) = 0
op = AffineFEOperator(a, l, U, V0)
uh = solve(op)
writevtk(Ω, "u_0_for_transcient", cellfields=["uh" => uh, "epsi" => ε(uh), "sigma" => σ ∘ ε(uh)])
uh
end
function solve_transcient(model, U0; order=1)
reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
V = TestFESpace(model, reffe;
conformity=:H1,
dirichlet_tags=["fixed"],
dirichlet_masks=[(true, true, true)]
)
g1_fixed(x, t) = VectorValue(0.0, 0.0, 0.0)
g1_fixed(t) = g1_fixed(0, 0)
U = TransientTrialFESpace(V, g1_fixed)
degree = 2 * order
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)
# try 5
m(t, u, v) = ∫(density * outer(v, u))dΩ
c(t, u, v) = ∫(0.0 * u ⋅ v)dΩ
a(t, u, v) = ∫(ε(v) ⊙ (σ ∘ ε(u)))dΩ
b(t, u, v) = ∫(0.0 * u ⋅ v)dΩ
b(t, v) = b(t, 0, v)
op = TransientAffineFEOperator(m, c, a, b, U, V)
linear_solver = LUSolver()
ode_solver = ThetaMethod(linear_solver, Δt, θ)
u_transcient = solve(ode_solver, op, U0, time_start, time_end)
# cut out saving to .vtu file to help with debugging.
function run_solver_t(uₕₜ)
for (uₕ, t) in uₕₜ
println(t)
end
end
run_solver_t(u_transcient)
end
model = make_model()
U0 = make_U0_solution(model)
U_trans = solve_transcient(model, U0)
Just as a quick reference to help out the conversation, here is the weak form on equation 4. https://opencfs.gitlab.io/userdocu/PDEExplanations/Singlefield/MechanicPDE/#governing-equations
OK. I changed line 90.
ode_solver = Newmark(linear_solver,Δt,0.5,0.25) # ode_solver = ThetaMethod(linear_solver, Δt, θ)
Now, I get
ERROR: This function belongs to an interface definition and cannot be used.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] macro expansion
@ C:\Users\garla\.julia\packages\Gridap\TyUsh\src\Helpers\Macros.jl:9 [inlined]
[3] solve_step!(uF::Vector{Float64}, solver::Newmark, op::Gridap.ODEs.TransientFETools.ODEOpFromFEOp{Gridap.ODEs.ODETools.Affine}, u0::Vector{Float64}, t0::Float64, cache::Nothing)
@ Gridap.ODEs.ODETools C:\Users\garla\.julia\packages\Gridap\TyUsh\src\ODEs\ODETools\ODESolvers.jl:16
[4] solve_step!(uF::Vector{Float64}, solver::Newmark, op::Gridap.ODEs.TransientFETools.ODEOpFromFEOp{Gridap.ODEs.ODETools.Affine}, u0::Vector{Float64}, t0::Float64)
@ Gridap.ODEs.ODETools C:\Users\garla\.julia\packages\Gridap\TyUsh\src\ODEs\ODETools\ODESolvers.jl:27
[5] iterate(sol::Gridap.ODEs.ODETools.GenericODESolution{Vector{Float64}})
@ Gridap.ODEs.ODETools C:\Users\garla\.julia\packages\Gridap\TyUsh\src\ODEs\ODETools\ODESolutions.jl:47
[6] iterate(sol::Gridap.ODEs.TransientFETools.TransientFESolution)
@ Gridap.ODEs.TransientFETools C:\Users\garla\.julia\packages\Gridap\TyUsh\src\ODEs\TransientFETools\TransientFESolutions.jl:65
[7] (::var"#run_solver_t#40")(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution)
@ Main C:\Users\garla\git\julia_playground\FEA\elastic_transcient_demo.jl:96
[8] solve_transcient(model::CartesianDiscreteModel{3, Float64, typeof(identity)}, U0::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}; order::Int64)
@ Main C:\Users\garla\git\julia_playground\FEA\elastic_transcient_demo.jl:101
[9] solve_transcient(model::CartesianDiscreteModel{3, Float64, typeof(identity)}, U0::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}})
@ Main C:\Users\garla\git\julia_playground\FEA\elastic_transcient_demo.jl:65
[10] top-level scope
@ C:\Users\garla\git\julia_playground\FEA\elastic_transcient_demo.jl:106
get_matrix
function because my operator is not linear. Any suggestion on how I can implement this? Thanks in advance!
I changed the inputs to the solver to have the U0, V0, A0.
ode_solver = Newmark(linear_solver,Δt,0.5,0.25) # ode_solver = ThetaMethod(linear_solver, Δt, θ)
V0 = interpolate_everywhere( VectorValue(0.0,0.0,0.0),U(0.0))
A0 = interpolate_everywhere( VectorValue(0.0,0.0,0.0),U(0.0))
u_transcient = solve(ode_solver, op, U0, V0, A0, time_start, time_end)
Now, I get this error.
ERROR: MethodError: no method matching Gridap.ODEs.TransientFETools.TransientFESolution(::Newmark, ::Gridap.ODEs.TransientFETools.TransientFEOperatorFromWeakForm{Gridap.ODEs.ODETools.Affine}, ::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, ::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, ::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, ::Float64, ::Float64)
Closest candidates are:
Gridap.ODEs.TransientFETools.TransientFESolution(::Gridap.ODEs.ODETools.ODESolver, ::TransientFEOperator, ::Any, ::Real, ::Real) at C:\Users\garla\.julia\packages\Gridap\TyUsh\src\ODEs\TransientFETools\TransientFESolutions.jl:13
Gridap.ODEs.TransientFETools.TransientFESolution(::Gridap.ODEs.ODETools.ODESolver, ::TransientFEOperator, ::Tuple, ::Real, ::Real) at C:\Users\garla\.julia\packages\Gridap\TyUsh\src\ODEs\TransientFETools\TransientFESolutions.jl:24
Gridap.ODEs.TransientFETools.TransientFESolution(::Any, ::Any) at C:\Users\garla\.julia\packages\Gridap\TyUsh\src\ODEs\TransientFETools\TransientFESolutions.jl:8
...
Stacktrace:
[1] solve(solver::Newmark, op::Gridap.ODEs.TransientFETools.TransientFEOperatorFromWeakForm{Gridap.ODEs.ODETools.Affine}, u0::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, v0::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, a0::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, t0::Float64, tf::Float64)
@ Gridap.ODEs.TransientFETools C:\Users\garla\.julia\packages\Gridap\TyUsh\src\ODEs\TransientFETools\TransientFESolutions.jl:51
[2] solve_transcient(model::CartesianDiscreteModel{3, Float64, typeof(identity)}, U0::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}; order::Int64)
@ Main C:\Users\garla\git\julia_playground\FEA\elastic_transcient_demo.jl:97
[3] solve_transcient(model::CartesianDiscreteModel{3, Float64, typeof(identity)}, U0::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}})
@ Main C:\Users\garla\git\julia_playground\FEA\elastic_transcient_demo.jl:65
[4] top-level scope
@ REPL[12]:1
All code with the updates is here https://gist.github.com/garland3/aba405ac597d3e3a3cca64d775516dc1.
transientfeoperator()
function. The error is at line::44 using Gridap
using Gridap.ODEs
using GridapDistributed
using PartitionedArrays
using Test
using BenchmarkTools
# using MPI
# GridapDistributed.PRange
function main(parts)
θ = 0.2
u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
u(t::Real) = x -> u(x,t)
f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)
domain = (0,1,0,1)
partition = (10,10)
model = CartesianDiscreteModel(parts,domain,partition)
order = 2
reffe = ReferenceFE(lagrangian,Float64,order)
V0 = FESpace(
model,
reffe,
conformity=:H1,
dirichlet_tags="boundary"
)
U = TransientTrialFESpace(V0,u)
Ω = Triangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)
#
m(u,v) = ∫(u*v)dΩ
a(u,v) = ∫(∇(v)⋅∇(u))dΩ
b(t,v) = ∫(v*f(t))dΩ
res(t,u,v) = a(u,v) + m(∂t(u),v) - b(t,v)
jac(t,u,du,v) = a(du,v)
jac_t(t,u,dut,v) = m(dut,v)
op = TransientFEOperator(res,U,V0;order = 1)
op_constant = TransientConstantMatrixFEOperator(m,a,b,U,V0)
t0 = 0.0
tF = 1.0
dt = 0.1
U0 = U(0.0)
uh0 = interpolate_everywhere(u(0.0),U0)
ls = LUSolver()
ode_solver = ThetaMethod(ls,dt,θ)
sol_t = solve(ode_solver,op,uh0,t0,tF)
sol_t_const = solve(ode_solver,op_constant,uh0,t0,tF)
# l2(w) = w*w
# tol = 1.0e-6
# for (uh_tn, tn) in sol_t
# e = u(tn) - uh_tn
# el2 = sqrt(sum( ∫(l2(e))dΩ ))
# @test el2 < tol
# end
# for (uh_tn, tn) in sol_t_const
# e = u(tn) - uh_tn
# el2 = sqrt(sum( ∫(l2(e))dΩ ))
# @test el2 < tol
# end
# createpvd("test-transient-heat-equation") do pvd
for (uh,t) in sol_t
writevtk(Ω,joinpath(@__DIR__,"..","mpi_test","test-transient-heat-equation_$t"),cellfields=["u"=>uh])
# println(evaluate(uh.fields.part,Point(0.5,0.5)))
# println(fieldnames(typeof(Ω.model.models.part.grid)))
end
# end
end
# end #module
partition = (2,2)
# @btime prun(main, mpi, partition)
prun(main, mpi, partition)
# MPI.Finalize()
# define integration measures
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)
Γ = BoundaryTriangulation(model)
dΓ = Measure(Γ, degree)
n = get_normal_vector(Γ)
# parameters for domain transformation and closure
y = [1.0, 2.0, 3.0]
J = 3
T(x) = T_fun(x, y, J)
transpDTinv(x) = TensorValue(DTinv_fun(x, y, J)')
# transpDTinvh = interpolate_everywhere(transpDTinv, V)
JT(x) = JT_fun(x, y, J)
# define variational problems
f(x) = κ0
g(x) = -κ0*1.0im
a(u, v) = ∫((transpDTinv ⋅ ∇(u)) ⋅ (transpDTinv ⋅ ∇(v)) * JT) * dΩ - ∫( κ0^2 * u * v * JT) * dΩ + ∫(im * κ0 * u * v * JT*norm(transpDTinv ⋅ n)) * dΓ
L(v) = ∫(f * v * JT) * dΩ - ∫(g * v * JT*norm(transpDTinv ⋅ n)) * dΓ
MultiFieldFESpaces
, but I am not sure if that is what I would use? Does anyone know of a good example like this, either in a test case or pull request, etc? I will continue to go through the tutorials, but perhaps I am looking in the wrong places.
Hello, I have been looking into the linear elasticity an-isotropic model tutorial, https://gridap.github.io/Tutorials/dev/pages/t003_elasticity/.
In my problem i have a weak form like ∫ (Young_module ∘ (E1, E2, tags) ⊙ ( Δ(v) ⋅ Δ(u) ))*dΩ
, where E1 and E2 are the two different youngs_module parameters, when I define the young's modulus function like below:
function Young_module(E1,E2,tags)
if tag == material1
return E1
elseif tag == material2
return E2
end
It shows me an error. I would like to know if this is the correct way to define the Young's modulus function that can take into account the spatial variation over the domain. if required I can also reproduce an MWE for convenience.
GridapDistributed.jl
, that has produced the same error that i was talking about :- It is not possible to perform the operation "inner" on the given cell fields.
Any help in this regard will be highly appreciated using Gridap.Geometry
using GridapDistributed
using PartitionedArrays
using Gridap
using GridapGmsh
function main(parts)
model = GmshDiscreteModel(parts,"/path/to/solid.msh")
order = 1
reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)
V0 = TestFESpace(model,reffe; conformity=:H1,dirichlet_tags=["surface_1","surface_2"], dirichlet_masks=[(true,false,false), (true,true,true)])
g1(x) = VectorValue(0.005,0.0,0.0)
g2(x) = VectorValue(0.0,0.0,0.0)
U = TrialFESpace(V0,[g1,g2])
degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)
labels = get_face_labeling(model)
labels = labels.labels.part
dimension = 3
tags = get_face_tag(labels,dimension)
alu_tag = get_tag_from_name(labels,"material_1")
function lame_parameters(E,ν)
λ = (E*ν)/((1+ν)*(1-2*ν))
μ = E/(2*(1+ν))
(λ, μ)
end
E_alu = 70.0e9
ν_alu = 0.33
(λ_alu,μ_alu) = lame_parameters(E_alu,ν_alu)
E_steel = 200.0e9
ν_steel = 0.33
(λ_steel,μ_steel) = lame_parameters(E_steel,ν_steel)
function σ_bimat(ε,tag)
if tag == alu_tag
return λ_alu*tr(ε)*one(ε) + 2*μ_alu*ε
else
return λ_steel*tr(ε)*one(ε) + 2*μ_steel*ε
end
end
a(u,v) = ∫( ε(v) ⊙ (σ_bimat∘(ε(u),tags)) )*dΩ
l(v) = 0
op = AffineFEOperator(a,l,U,V0)
uh = solve(op)
writevtk(Ω,"demo_bimat",cellfields=["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ_bimat∘(ε(uh),tags)])
end
partition = (1,3,2)
prun(main, mpi, partition)
Hi all, I am having problems precompiling GridapMakie after updating to Julia 1.8.1. Any guidance? Clue? This is the error I get:
``(base) reula@boga Tarea_5 % julia
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _
| |
| | || | | | (| | | Version 1.8.1 (2022-09-06)
/ |__'|||_'| | Official https://julialang.org/ release
|__/ |
julia> import Pkg; Pkg.precompile()
Precompiling project...
✗ GridapMakie
0 dependencies successfully precompiled in 71 seconds. 465 already precompiled.
ERROR: The following 1 direct dependency failed to precompile:
GridapMakie [41f30b06-6382-4b60-a5f7-79d86b35bf5d]
Failed to precompile GridapMakie [41f30b06-6382-4b60-a5f7-79d86b35bf5d] to /Users/reula/.julia/compiled/v1.8/GridapMakie/jl_82qz3m.
ERROR: LoadError: UndefVarError: point_iterator not defined
Stacktrace:
[1] getproperty(x::Module, f::Symbol)
@ Base ./Base.jl:31
[2] top-level scope
@ ~/.julia/packages/GridapMakie/oS6vj/src/recipes.jl:158
[3] include(mod::Module, _path::String)
@ Base ./Base.jl:419
[4] include(x::String)
@ GridapMakie ~/.julia/packages/GridapMakie/oS6vj/src/GridapMakie.jl:1
[5] top-level scope
@ ~/.julia/packages/GridapMakie/oS6vj/src/GridapMakie.jl:17
[6] include
@ ./Base.jl:419 [inlined]
[7] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::Nothing)
@ Base ./loading.jl:1554
[8] top-level scope
@ stdin:1
in expression starting at /Users/reula/.julia/packages/GridapMakie/oS6vj/src/recipes.jl:158
in expression starting at /Users/reula/.julia/packages/GridapMakie/oS6vj/src/GridapMakie.jl:1
in expression starting at stdin:1
Stacktrace:
[1] pkgerror(msg::String)
@ Pkg.Types /Applications/Julia-1.8.1.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/Types.jl:67
[2] precompile(ctx::Pkg.Types.Context, pkgs::Vector{String}; internal_call::Bool, strict::Bool, warn_loaded::Bool, already_instantiated::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Pkg.API /Applications/Julia-1.8.1.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/API.jl:1427
[3] precompile
@ /Applications/Julia-1.8.1.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/API.jl:1058 [inlined]
[4] #precompile#225
@ /Applications/Julia-1.8.1.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/API.jl:1057 [inlined]
[5] precompile (repeats 2 times)
@ /Applications/Julia-1.8.1.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/API.jl:1057 [inlined]
[6] top-level scope
@ REPL[1]:1
julia>
```
Hi all, I used to plot the real part of the first component of a complex vector using
fig, ax, plt = plot(Ω, ( u-> real(u[1])) ∘ uₕ)
But now I am getting an error which I guess it is related to dimensions:
Error: error handling request
│ exception = (ErrorException("Metadata array needs to have same length as data.\n Found 2400 data items, and 1600 metadata items”)
include("ConstantFESpaces.jl")
is missing in the FESpaces.jl file.
interpolate_everywhere
in each non-linear operation (*, /, and others). I would like to know whether there are others ways. Many thanks in advance.
-pc_type fieldsplit -pc_fieldsplit_type schur
as a preconditioner for a multifield problem using GridapPETSc? In the petsc options, for the DarcyDriver example, I am using the following preconditioner -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition a11 \
-fieldsplit_0_pc_type lu \
-fieldsplit_1_ksp_rtol 1.e-9 -fieldsplit_1_pc_type lu
It does not recognize that there are 2 fields. Thanks in advance
CartesianDiscreteModel
to do a 3d mesh, and I need to give diferent boundary conditions on each face of the box. But im confuse to how do it, i dont understand what part of the box are the labels "tag_01"- "tag_27" to build the diferent faces. II need to understand the labels. If someone could help me i will apresiate it
uh = solve(op) #solving some problem
l2(v)= ∫(exp(uh)*v)dΩ
op2 = AffineFEOperator(a,l2,U,V)
uh2 = solve(op)
disp_max = 43
disp_inc = 1
nsteps = ceil(Int,abs(disp_max)/disp_inc)
x0 = zeros(Float64,num_free_dofs(V0_Disp))
cache = nothing
for step in 1:nsteps
disp_x = step * disp_max / nsteps
x0, cache,uh = run(x0,disp_x,step,nsteps,cache)
Force = sum(∫( n_Γ_Load ⋅ (S∘(∇(uh))) ) *dΓ_Load)
push!(Load, Force[2])
push!(Displacement, disp_x)
end