Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Activity
  • 03:08
    codecov-commenter commented #838
  • 02:05
    JordiManyer synchronize #838
  • 02:05

    JordiManyer on refined-discrete-models

    Added get_refined_polytopes (compare)

  • Dec 05 22:57
    hgpeterson closed #850
  • Dec 05 22:57
    hgpeterson commented #850
  • Dec 05 16:57
    codecov-commenter commented #852
  • Dec 05 16:26
    JordiManyer synchronize #838
  • Dec 05 16:26

    JordiManyer on refined-discrete-models

    Compressed is_change_possible a… (compare)

  • Dec 05 15:51
    pmartorell edited #852
  • Dec 05 15:50
    pmartorell assigned #852
  • Dec 05 15:50
    pmartorell review_requested #852
  • Dec 05 15:50
    pmartorell opened #852
  • Dec 05 15:29
    JordiManyer synchronize #838
  • Dec 05 15:29

    JordiManyer on refined-discrete-models

    Small change to a docstring (compare)

  • Dec 05 15:26
    JordiManyer synchronize #838
  • Dec 05 15:26

    JordiManyer on refined-discrete-models

    Compacted is_related code (compare)

  • Dec 05 15:18
    JordiManyer synchronize #838
  • Dec 05 15:18

    JordiManyer on refined-discrete-models

    Added check in AdaptedTriangula… (compare)

  • Dec 05 14:41
    JordiManyer synchronize #838
  • Dec 05 14:40

    JordiManyer on refined-discrete-models

    Added a new adapt method for fu… (compare)

Paul E. Mendez
@Paulms
Hello, I need to add to my weak formulation a parameter that depends on a norm of a Tensor variable (i.e. a non-linear viscosity for a viscoplastic fluid). I tried to use a ternary operator like this: 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
? Thanks.
2 replies
Devasmit Dutta
@DevasmitDutta
hello, I would like to know if this is the correct way, of evaluating a ::DistributedCellfield at any point on the mesh. For example, println(evaluate(uh.fields.part,get_cell_coordinates(Ω.trians.part)[1][2])), returns an error
Error: │ exception = │ AssertionError: A check failed
I am just wondering if its due to any bug, because on a CartesianDiscreteModel() this works fine, but when I try to evaluate on a un-structuredDiscretemodel() it shows A check failed error.
2 replies
Stauber López Daniela
@DaniStauber
Hi everybody! Thanks for the dynamics of this foro, helps it a lot!! I wanna use the tool of FESpaceWithLinearConstraints, but im a little lost about how use it. If anyone has seen implementations or found any documentation that can shed light on this, I'd really appreciate it.
3 replies
Anthony Garland
@AnthonyPGarland_twitter

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)
Oriol Colomés
@oriolcg
Hi @AnthonyPGarland_twitter what error do you get? Just looking at the code I see that you use Thetamethod, which doesn’t work for 2nd order ODEs. Try using Newmark.
Anthony Garland
@AnthonyPGarland_twitter

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
Oriol Colomés
@oriolcg
Hi @AnthonyPGarland_twitter , Newmark method requires initial conditions for displacement, velocity and acceleration. You are only sending displacement.
carlodev
@carlodev
Hi, I am solving a non-linear, time dependent, multifield problem using PETSc as solver, in particular using gamg as preconditioner. As mentioned in petsc documentation you have to set MatSetBlockSize(), as in this test https://github.com/gridap/GridapPETSc.jl/blob/master/test/sequential/ElasticityDriver.jl . However, I cannot use get_matrix function because my operator is not linear. Any suggestion on how I can implement this? Thanks in advance!
4 replies
philippweder
@philippweder
Hello!
I am quite new to Gridap.jl and I have a question regarding the general syntax rules of weak formulations. More precisely, I want to write a weak form, where I map my problem back to a reference domain. As a consequence, the Jacobian matrices appear in my integrands. How can I include them without getting errors?
Anthony Garland
@AnthonyPGarland_twitter

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.

5 replies
Devasmit Dutta
@DevasmitDutta
Hello, I am attaching a below MWE, where I am trying to implement auto-diff in GridapDistributed code. It will be really good if there is any-way, to utilise auto-diff, using GridapDistributed code here. As below code, shows error when trying to do that, but works fine, when I am explicitly passing both residual and jacobian into the 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()
3 replies
philippweder
@philippweder
How can I implement the expression ||DT^{-1} n|| in a variational formulation, where DT is the explicit implementation of the Jacobian matrix of a function and n is the outward normal vector on the boundary? The following code throws an error, which I don't manage to decipher:
# 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)) *
Krishna Bhogaonker
@00krishna
Hey folks, I am totally new to Gridap, so just getting used to the package. But i was wondering if there is any feature comparison between Gridap and say Fenics or other packages. I don't want to migrate a bunch of code until I know that Gridap has everything that I would need. Is there like a table or such? Or even easier, are there any big feature gaps between Gridap and packages like Fenics or FreeFem, etc?
7 replies
CarlosContrerasQ12
@CarlosContrerasQ12
Hi, I am trying to get an array of the my solution evaluated at my nodes, but I think using evaluate(u,Point(node)) is not very efficient, ¿are the values of my solution stored somewhere in my SingleFieldFEFunction() object? I tried the solution proposed here gridap/Gridap.jl#345 , but didn't work as it seems the function get_cell_values(u) is not implemented
10 replies
Krishna Bhogaonker
@00krishna
hey folks, I have been going through the tutorials and they are very helpful. I was looking at a specific problem, and I am not sure if there is a tutorial that addresses it? So say I have an elastic wave from an earthquake, and I want to see how this wave propagates through different layers of rock. So the rock strata can have different permeability constants depending on density, etc. So I imagine I would look at this like having 2 different meshes--one for each strata, and then the two meshes are joined at the boundary point. In the tutorials I found some examples of hyperbolic problems and these 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.
32 replies
Roman
@roman-schaerer
Hi, I would like to implement some stabilization method for which I need to acquire the local element size. I'm not yet familiar with the more lower-level API calls in Gridap that might be required to compute e.g. the cell-wise element size, e.g. the longest or shortest edge length. What would be a possible approach to this task?
6 replies
Antoine Marteau
@Antoinemarteau
Hi guys :) Has someone implemented a gauge condition for a (unknown) vector potential in Nedelec space ?
I'm trying the co-tree gauge, with tree build by gmsh, is a dirichlet constraint appropriate to set inner tree edges DOF to 0 ? Should the dirichlet function return a scalar 0 or a VectorValue'd 0. ? Thanks !
2 replies
Devasmit Dutta
@DevasmitDutta

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.

3 replies
Oriol Colomés
@oriolcg
What is the error that you get? Is tag a cellfield?
5 replies
Devasmit Dutta
@DevasmitDutta
Hello @oriolcg , below is the MWE on the same linear-elasticity bi-material tutorial code that has been parallelised using 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*μ_aluelse
              return λ_steel*tr(ε)*one(ε) + 2*μ_steelend
      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)
philippweder
@philippweder
Are complex-valued Nedelec elements not available in Gridap? I would like to solve a time-harmonic Maxwell problem. Is there some workaround?
13 replies
55113110039
@55113110039
Hi, I am new to Gridap. Can I define tags on plain a CartisianDiscreteModel without the help of a model.json or model.gsh file?
2 replies
Sounak Das
@Sounakdas013
how to impose plane stress or plane strain on the hyperelasticity code given in the tutorial?
1 reply
Sounak Das
@Sounakdas013
Can anyone confirm if the hyperelasticity tutorial code is entirely correct or not?
2 replies
philippweder
@philippweder
In code validation, how to correctly impose Neumann boundary conditions for a known solution?
8 replies
Eric Neiva
@eneiva_gitlab
:loudspeaker: :loudspeaker: :loudspeaker: Hello everybody! @fverdugo and I are going to teach a Gridap tutorial by invitation of Groupe Calcul in :pushpin: Paris (INRIA Saclay) next :calendar: December 1! Registration is open :point_right: https://indico.mathrice.fr/event/365/
Looking forward to seeing you there!
9 replies
Oscar Reula
@reula

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>
```

3 replies
Z J Wegert
@zjwegert
Hi all, I'm trying to compute N integrals which have dependence on the shape function for each node (for i = 1:N). E.g., vᵢ = ∫ f(v)ϕᵢ dΩ for i = 1:N. Is there a nice way of doing this kind of calculation in Gridap? Cheers, Zach
ziolai
@ziolai
I am trying to define a convection-diffusion operator. Using v(x) = VectorValue( 1.0, 0.0 )
a(u,v) = ∫( ∇(v)⋅∇(u) + ∇(u)⋅v)*dΩ fails. Pls. advice.
2 replies
philippweder
@philippweder
I am trying to use GridapPETSc.jl in a serial way as indicated here. However, my matrices are complex-valued. How can I tell PETSc to use complex numbers?
17 replies
carlodev
@carlodev
Hello everyone, I have developed a package for automatizing and optimising the procedure for creating structured airfoil meshes using Gmsh and creating physical boundaries. It is far from perfect, but it may help some of you. if you want to give it a look https://github.com/carlodev/GmshAirfoil ; feel free to contribute
philippweder
@philippweder
Maybe a more general question: I am solving a 3D time-harmonic Maxwell problem, and I would like to improve the performance (on a single machine). What can I do to achieve that?
I have added a SparseMatrixAssembler. Does this really change anything? Which linear solver should I use?
ziolai
@ziolai
Helmholtz or full Maxwell?
philippweder
@philippweder
I am actually doing both, but I am only struggling with the full Maxwell problem
Paul E. Mendez
@Paulms
Hello, After solving a problem I get a vector variable u (i.e. u = solve(AffineFEOperator(a,L,X,Y))), How can I compute the L infinite norm of the divergence of u? I can get the dof values for u, but not sure how to get values of div(u).
Oscar Reula
@reula

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”)
philippweder
@philippweder
Are there ways to effectively reduce the memory footprint of Gridap simulations? My server crashed multiple times because Gridap seems to allocate tons of memory...
More precisely, although I am using a sparse assembler, Gridap seems to instantiate a dense matrix before solving the linear system. Is there a way to prevent that?
11 replies
Paul E. Mendez
@Paulms
Hello, the changelog for the last version, mention that a constant space has been added, however it seems that it is not being exported. I tried to follow the test, but i got "ERROR: UndefVarError: ConstantFESpace not defined". I wonder if a include("ConstantFESpaces.jl") is missing in the FESpaces.jl file.
Sounak Das
@Sounakdas013
Is finite deformation problems very much dependent on mesh size and the displacement increment step applied on the boundary- I am getting some sharp spike in the load-displacement curve. Can anyone help me with that ?
hyun825
@hyun825
Hello, I am working in a convection-difusion-reaction PDE coupled to some first order non-linear ODEs through the reaction terms of PDE. I have implemented a Newton-Raphson for the ODEs which computes the variables in the FE space in order to solve the time-dependent PDE. I have tried in different ways, and the way I have found which consumes less compilation time and less computation time was by using interpolate_everywhere in each non-linear operation (*, /, and others). I would like to know whether there are others ways. Many thanks in advance.
Oscar Reula
@reula
Hi, I am trying an eigenvector problem where I have a surface term: ∫( (u×n_n)⋅v )*dΓ but I get an error related to conversion to or from Complex64. I saw this was also in an old question, but I didn’t see the solution on the thread. Please help!
5 replies
carlodev
@carlodev
Has anyone used -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
4 replies
guin0x
@guin0x
does anyone have an example of a damage model for multi-materials? (following something like this)
Stauber López Daniela
@DaniStauber
Hi there! I want to use the 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
Francesc Verdugo
@fverdugo
@daniela write the cartesian model in vtu format and open the resulting files with paraview. There you can inspect the labels
1 reply
Nanna Berre
@nannaberre
Hi! I want to use the exponential of a previous solution. Is there any way for the current functionality to do the following:

uh = solve(op) #solving some problem

l2(v)= ∫(exp(uh)*v)dΩ
op2 = AffineFEOperator(a,l2,U,V)
uh2 = solve(op)

1 reply
Sounak Das
@Sounakdas013
I want to calculate load at a boundary for a finite deformation problem. sum(∫( n_Load ⋅ (S∘(∇(uh))) ) *dΓ_Load. Will this equation give me the correct value of load as n_Load is the normal vector at the boundary in the reference configuration and S is the second piola kirchoff tensor. Using software like comsol I am getting same value in literature but using hyperelasticity code the value in load-displacement curve it is differing a lot. If anyone can clarify , it will be helpful.
3 replies
Sounak Das
@Sounakdas013
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
can anyone confirm whether the load calculation is correct or not or if any modifications in the equation are needed ?