Gridap.jl/community

Santiago Badia

@santiagobadia

In this tutorial we have tried to explain the key design concepts in Gridap, e.g.., lazy maps and arrays, etc/

Take a look and let us know what you think

Bhavesh Shrimali

@bhaveshshrimali

Thanks @santiagobadia . I was about to point to that (but wasn't sure if it is the best first step) especially since I find myself taking a look at the dev tutorial very often. This in addition to peeking into the repository and the branches, and experimenting in REPL has cleared up a lot of issues in the past for me.

Erik Schnetter

@eschnett

thanks! that's very useful.

i see that a cartesian model has two important fields, called grid and grid_topology. the grid seems to be defined in an abstract reference space, and the grid_topology seems to contain the physical geometry of the grid. is that correct? if so, the name grid_topology seems to be a misnomer. (no criticism here.)

the abstract grid is also called "triangulation", i.e. a triangulation can consist of cubes, and it doesn't have coordinates attached.

Francesc Verdugo

@fverdugo

Grid has information abou nodes (i.e the geometrical map) and grid topology has info about vertices, edges, faces, cells and their incidence relation. Grid and triangulation are overlaped and will be possibly merged in the future.

Paul E. Mendez

@Paulms

Hi! I am trying to solve a nonlinear problem (p-Laplacian in mixed formulation). In the p-Laplacian tutorial you can call norm∘grad(u), but replacing grad(u) for a Tensorial Variable doesnt seem to work. Adding the method: norm(u::MultiValue{Tuple{D1,D2}}) where {D1,D2} = sqrt(inner(u,u)) solves part of the problem (I can call norm∘T with T being tensor), but (norm∘T)^(p-2), still errors. Any clue of what other method I would need to add for this to work. Thanks!

Zach Wegert

@Omega-xyZac

I may be wrong but I believe norm is already defined for 2-tensors. To sort out the second part of your question, you should define a function: ‘p_power(t)=t^(p-2)’ then call ‘p_power ∘ (norm ∘ T)’

For some reason I can’t get code to work on my phone. Apologies

That might do the trick

Paul E. Mendez

@Paulms

Yeah, thanks!!!, that worked. Though, It's not clear for me when to use p_power(norm(T)) and when p_power∘(norm(T)) . In other hand, If I define the space as TestFESpace(model,ReferenceFE(lagrangian,TensorValue{2,2,Float64,4},order),conformity=:L2) the norm is not defined. it works adding the method of my previous message.

Zach Wegert

@Omega-xyZac

Great! I believe you should be using p_power∘(norm(T)) when you are applying p_power to an array of nodal data - @fverdugo would probably answer this question better. It appears that this is not yet defined in Operations.jl for MultiValue{Tuple{D1,D2}} but is for MultiValue{Tuple{D}}. May be worth adding.

Paul E. Mendez

@Paulms

Ok, should I open an Issue?

Zach Wegert

@Omega-xyZac

Sure, open an issue.

Paul E. Mendez

@Paulms

Ok, I added the issue: gridap/Gridap.jl#536

Romain Veltz

@rveltz

Hi, is there an easy way to impose periodic boundary conditions (basically a circle)?

Zach Wegert

@Omega-xyZac

Hi Romain, I believe PBCs are only implemented for Cartesian mesh’s. I think there is an open issue for non-Cartesian mesh’s

Romain Veltz

@rveltz

I am happy with cartesian mesh! My issue is how to specify it for example in V:

using Gridap
n = 40
domain = (0,1)
cells = (n,)
model = CartesianDiscreteModel(domain,cells)
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe,conformity=:H1,)

Romain Veltz

@rveltz

For example, this does not work

using Gridap

u(x) = sin(2*pi*(x[1]-0.25))*x[2]
f(x) = -Δ(u)(x)

domain = (0,1,0,2)
cells = (10,20)
model = CartesianDiscreteModel(domain,cells;isperiodic=(true,false))

Francesc Verdugo

@fverdugo

I believe PBCs are only implemented for Cartesian mesh’s.

GridapGmsh already supports periodic BCs!

using Gridap
u(x) = sin(2pi(x[1]-0.25))*x[2]
f(x) = -Δ(u)(x)
domain = (0,1,0,2)
cells = (10,20)
model = CartesianDiscreteModel(domain,cells;isperiodic=(true,false))

This code works just fine in my repl

make sure that you are using a recent version

Romain Veltz

@rveltz

Are you on 1.5? For

(@v1.6) pkg> st Gridap
      Status `~/.julia/environments/v1.6/Project.toml`
  [56d4f2e9] Gridap v0.15.1 `https://github.com/gridap/Gridap.jl.git#master`

The above example returns

julia> model = CartesianDiscreteModel(domain,cells;isperiodic=(true,false))
ERROR: BoundsError: attempt to access 1-element Vector{Int32} at index [0]
Stacktrace:
  [1] getindex
    @ ./array.jl:801 [inlined]
  [2] _face_to_cells_count!
    @ ~/.julia/packages/Gridap/zovI1/src/Geometry/GridTopologies.jl:889 [inlined]
  [3] _face_to_cells(cell_to_faces_data::Vector{Int64}, cell_to_faces_ptrs::Vector{Int32}, nfaces::Int64)
    @ Gridap.Geometry ~/.julia/packages/Gridap/zovI1/src/Geometry/GridTopologies.jl:867
  [4] generate_cells_around
    @ ~/.julia/packages/Gridap/zovI1/src/Geometry/GridTopologies.jl:564 [inlined]
  [5] Gridap.Geometry.UnstructuredGridTopology(vertex_coordinates::Vector{VectorValue{2, Float64}}, cell_vertices::Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}, cell_type::Vector{Int8}, polytopes::Vector{Gridap.ReferenceFEs.ExtrusionPolytope{2}}, orientation_style::Gridap.Geometry.Oriented)
    @ Gridap.Geometry ~/.julia/packages/Gridap/zovI1/src/Geometry/UnstructuredGridTopologies.jl:36
  [6] _generate_grid_topology_from_grid(grid::Gridap.Geometry.UnstructuredGrid{2, 2, Float64, Gridap.Geometry.Oriented}, cell_to_vertices::Gridap.Arrays.Table{Int64, Vector{Int64}, Vector{Int32}}, vertex_to_node::Vector{Int64})
    @ Gridap.Geometry ~/.julia/packages/Gridap/zovI1/src/Geometry/UnstructuredGridTopologies.jl:148
  [7] _cartesian_grid_topology_with_periodic_bcs
    @ ~/.julia/packages/Gridap/zovI1/src/Geometry/CartesianDiscreteModels.jl:413 [inlined]
  [8] CartesianDiscreteModel(desc::Gridap.Geometry.CartesianDescriptor{2, Float64, typeof(identity)})
    @ Gridap.Geometry ~/.julia/packages/Gridap/zovI1/src/Geometry/CartesianDiscreteModels.jl:20
  [9] #CartesianDiscreteModel#82
    @ ~/.julia/packages/Gridap/zovI1/src/Geometry/CartesianDiscreteModels.jl:69 [inlined]
 [10] top-level scope
    @ REPL[13]:1

Francesc Verdugo

@fverdugo

I am using Julia 1.5

Alberto F. Martín

@amartinhuertas_gitlab

Hi @rveltz. We could reproduce the issue that you are reporting with Julia v1.6.0-rc1 in GH actions CI. See https://github.com/gridap/Gridap.jl/runs/1875200999?check_suite_focus=true#step:6:143. We will look into the issue and solve it.

Dominik Büschgens

@bueschgens

Hi all, I am new to Gridap and trying to solve the transient heat equation in three dimensions for an arbitrary mesh. I found an example with GridapODEs Plugin. In the beginning of the example there is an equation given u(x,t), I guess this is the solution, but then I do not understand the sense of the example. Since there is no documentation it would be nice if someone can explain what is happening there. Thank you in advance.

oriolcg

@oriolcg

Welcome to Gridap @bueschgens. The test that you mention solves a problem with a manufactured analytical solution. That means that you enforce a forcing term such that the solution satisfy a given function. In that case u(x,t) is the analytical function (or exact solution) and the forcing term is computed using this function: f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x). The Finite Element solution at time tn, uh_tn, is obtained when you iterate over the result of the solve function:

sol_t = solve(solver,op,uh0,t0,tF)

l2(w) = w*w

tol = 1.0e-6
_t_n = t0

for (uh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  e = u(tn) - uh_tn
  el2 = sqrt(sum( ∫(l2(e))dΩ ))
end

Note that you can also use the analytical function to compute the numerical error.

Erik Schnetter

@eschnett

when i calculate the gradient of a SingleFieldFEFunction via ∇(fh), then the result is a GenericCellField. how do i find out what FE space this field uses?

is there a way to "force" the gradient to be an FEFunction?

i encounter this when i try to implement evaluate! for the gradient. i'm currently implementing evaluate!(cache,f::SingleFieldFEFunction,x::Point). is that the right signature?

Francesc Verdugo

@fverdugo

Hi @eschnett ! I would implement this function for all CellFields not only for SingleFieldFEFunction. I think that the API of CellField is rich enough to implement this

Erik Schnetter

@eschnett

how do i obtain the FESpace for a CellField? get_fe_space isn't implemented e.g. for a GenericCellField.

Francesc Verdugo

@fverdugo

From a CellField cf, you can get_triangulation(cf). I think this is sufficient to implement evaluation at a point x.

Do you need the FESpace for some reason?

Erik Schnetter

@eschnett

i need to evaluate the function within an element. for that i get the FE space, then the FE basis, then i can evaluate there.

i already found the cell where i need to evaluate, now i need to evaluate within that cell.

Francesc Verdugo

@fverdugo

already found the cell where i need to evaluate

Great! This is the difficult part

once you have this @eschnett you just need to to some think like:

function return_cache(f::CellField,x::Point)
  trian = get_triangulation(f)
  setup = _setup_kd_tree_search(trian)
  cell_f = get_array(f)
  c1 = array_cache(cell_f)
  f = testitem(cell_f)
  c2 = return_cache(f,x)
  setup, c1, c2, cell_f, f
end

function evaluate!(cache,f::CellField,x::Point)
  setup, c1, c2, cell_g, g = cache
  if f === g
    cell_f = cell_g
  else
    cell_f = get_array(f)
  end
  cell = _find_cell_kd_tree(setup,x)
  f = getindex!(c1,cell_f,cell)
  fx = evaluate!(c2,f,x)
  fx
end

_setup_kd_tree_search represents a function that set-ups the kd-tree searches. _find_cell_kd_tree does the search for a given point

From what you say, it is likely that you have already implemented these functions ore similar ones

Erik Schnetter

@eschnett

thanks. yes, i have almost everything in place. the missing piece is getting the topology. i believe i can get a topology from the model, and i see that a model is a triangulation. i can get the triangulation, but when i construct a topology from it, then this is always an UnstructuredGridTopology. you hint to call get_grid_topology doesn't work since it requires a model, which i don't have, i only have a triangulation.

is there a generic function to check whether a cell contains a point? i'm sure there is a function to evaluate a basis function on a cell, and that condition should fall out on the side.

finally, i am encountering

Function inverse_map is not implemented yet for objects of type Gridap.Fields.LinearCombinationField{Vector{VectorValue{2, Float64}}, Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{2, Float64}}}

how do i solve this?

Erik Schnetter

@eschnett

@fverdugo thanks for your assistance earlier today. your pointers and pull requests helped. i think things are working now!

Erik Schnetter

@eschnett

i have a CartesianDiscreteModel, and i want to apply different boundary conditions on different faces of a square. how do i approach this? is there a way to get tags for BoundaryTriangulation?

oriolcg

@oriolcg

You can check tests/GridapTests/SurfaceCouplingTests.jl

Where communities thrive

People

Repo info

Activity